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I. AN INTRODUCTION TO ASSET 


The Antisubmarine Warfare Systems Evaluation Tool (ASSET) is a campaign-level 
simulation which models open ocean scenarios involving submarines, maritime patrol 
aircraft (MPA), shore-based command and data-fusion centers and a wide variety of 
passive acoustic sensors. Active acoustic and non-acoustic sensors are only modeled 
using a simple area sensor model which has a fixed probability of detection and false 
alarm rate. The simulation scenario is specified by a user-supplied "architecture" which 
determines all facets of environment, command control, sensor interaction and platform 
maneuver. This structure is input to a desktop computer/workstation (the Apple 
Macintosh II in this version) through a series of user-friendly windows, each dealing with 
a specific topic. A particular configuration of sensors, data fusion centers, 
communication nodes and tactical platforms interact in a particular geographic location 
against an analogous enemy force structure. 

The scenario is repeated as a Monte Carlo simulation to produce statistically 
meaningful measures of effectiveness (MOE’s). Output data regarding the detection, 
localization and prosecution of enemy submarines can then provide a quantitative basis 
for decisions regarding ASW Master Plans, Top Level Warfare Requirements, and a 
variety of emerging technology assessments and system appraisals [Ref. 1:p. iv]. These 


results can also be useful in conducting quantitative experimentation with regional force 


levels, force compositions and commitment strategies; command, control, 
communications and intelligence (C*I) networks; implications of foreign technology 
advances; and fleet exercise planning. The current version of ASSET (version 1.0) limits 
the scope of these scenarios to open ocean search and prosecution of hostile submarines 
by cueing friendly submarines and MPA with wide-area sensors, and modeling the 
supporting C*I networks. The modular nature of the object-oriented structure of ASSET 
makes expanding the scope of the simulation possible [Ref. 2:p. 1-2]. As the 
construction and operating costs of naval platforms spiral upwards, a simulation tool of 
this kind will be increasingly important to intelligently manage resources vis-a-vis a 
rapidly evolving spectrum of possible threats. 
A. OBJECTIVES AND RATIONALE 
The purpose of this thesis is three-fold: 
® to provide a description of the target tracking algorithm used by ASSET and 
how it relates to: 
1. real-world tracking. 
2. the ASSET simulation as a whole. 
® to examine the mathematical development of the tracking algorithm. 
® to examine possible modifications to the existing Kalman filter tracking 


algorithm to better match the data input to it and increase computational 
efficiency. 


The use of full Kalman filter equations for elliptical areas of uncertainty 
(AOU’s) is not necessary given the circular AOU’s that ASSET’s sensor model 
produces. Making use of the fact that the computer knows the ground truth locations 
of the platforms, and the limited data input to the filter, it will be shown that it is 
possible to approximate appropriately distributed filtered state positions using 
simplified forms of the Kalman filter matrix computations. The situations when the 
tracker algorithm is used can also be limited based on the information content of the 
contact. Together, these modifications yield a significant reduction in the 
computational complexity of the tracking algorithm, which is the principle goal of this 
thesis. 

The automatic correlator tracker (ACT) implemented in ASSET is a stand-alone 
module based on the Ocean Surveillance Information System (OSIS) baseline upgrade 
single hypothesis, multiple target, Kalman filter-based, correlation and tracking 
algorithm [Ref. 3:p. 4]. This set of routines contains a complete representation of the 
OSIS baseline upgrade automatic correlator tracker’s (OBU-ACT) functions, with the 
exception of LINK, the utility that evaluates track sets which are potentially legs of a 
single target’s track; and EQUATE, which makes track associations for contacts pre- 
correlated to a particular platform based on acoustic or electromagnetic signature. 

The OBU-ACT package is designed to process contact information derived from 
position-only, bearing-only, and position and velocity reports. The covariance 
matrices associated with these reports can be interpreted as an elliptical representation 


of the error in position and velocity. The capability to process line of bearing 


contacts is not available in ASSET (1.0). The method used to model sensor systems 
in ASSET 1.0 approximates the elliptical errors with circular ones. This reduces the 
complexity of the covariance matrices considerably. Taking advantage of the large 
number of zero elements and repeated values in the matrices actually constructed by 
ASSET, the matrix equations involved in the tracking algorithm can be reduced to 
equivalent scalar ones by eliminating all the zero-multiplied reer 

Central to this algorithm is the Integrated Ornstein-Uhlenbeck (IOU) process 
which models target motion in the Maneuvering Target Statistical Tracker (MTST). 
The MTST uses this model, together with contact data corrupted by noise, to estimate 
the size and position of the Submarine Probability Area (SPA) which has an 86 
percent probability of containing the target’s true position. The major modification to 
the algorithm proposed is based on assuming that the IOU prediction position 
distribution, in like fashion to the contact position distribution, is centered on the 
target’s true position. This results in a distribution of SPA centers that is centered on 
the targets true position also. Thus, only the SPA variance is calculated, not its 
center position, which is drawn randomly from the estimated distribution about the 
true target position. Since contact data does not play a role in computing these 
variances, they may be calculated prior to the start of the simulation. A table of 
variances for the sensor AOU values defined for the scenario and an appropriate 
range of contact report interarrival times can be constructed and referred to as needed 


during simulation execution to minimize the time spent processing contact data. 


B. FUNDAMENTALS OF OPERATION. 


The basic structure of ASSET was adapted from COAST, the Common LISP 
Architectural Study Tool. Common LISP is an object-oriented programming language 
used for rapid prototyping and artificial intelligence work. A general overview of this 
structure is presented to facilitate understanding of how the COAST/ASSET system is 


organized and the critical role of the tracker to simulation operation. 


1. Objects and Object Interactions. 

Object-oriented programs combine data with associated procedures to form a 
hierarchical network of self-contained modules known as objects. These objects can 
then be invoked by each other according to the program methodology. Program 
objects are independent and can be modified without affecting other objects they 
communicate with. The data that is communicated will change, but the relationship 
between the objects does not necessarily have to. Data can be evaluated internally by 
an object without affecting the analogous data in a separate object. 

The principle of instancing allows a single object definition to create any 
number of structures within the program, all obeying the same definition. This 
allows a single definition of a submarine object, for example, to be given many 
different sets of parameters each representing a different class of submarine. Other 
objects may be designated to impart their characteristics to the new submarine object. 
This process of inheritance allows the properties of another object, an acoustic sensor 


object for example, to be created separately from the submarine object. The 


submarine object can then be designated to be a kind of acoustic sensor and it will 
inherit the properties of the acoustic sensor object in addition to its own. Another 
object such as a SURTASS towed-array surveillance ship object can also inherit the 
properties of the acoustic sensor object, with different parameters, without affecting 
the parameters of the acoustic sensor inherited by the submarine object. 

Complex classes of platforms can be formed by multiple inheritance of the 
properties of basic building block objects. These are then used to simulate any 
number of platforms of that class, each operating independently within the simulation. 
The computer memory available to store each of these instances of the object becomes 
the limiting factor to the complexity of the scenario to be examined. The organization 
of objects in ASSET consist of several major groups: 

@ Graphical environment objects which implement the standard Macintosh 
graphical user interface to provide input/output by means of a mouse based 
"point and click" metaphor featuring pull-down menus,dialog boxes with labeled 
data regions and radio buttons, and windowed presentation of multiple 
information sources simultaneously. 

@ Region management objects which construct and manage distinct regions of 
varying environmental properties, command responsibility or platform patrol 


assignment. 


@ Event management objects which queue all simulation events in time order 
sequence and parcel them out to the appropriate resolution objects. 


@ Command objects which perform resource allocation, assigning available assets 
to contact cues based on either time to station (MPA) or area of uncertainty size 
(submarines). 


@ Automatic Correlator Tracker (ACT) objects which manage the grouping and 
processing of contact reports into target tracks creating a tactical picture 
consisting of combinations of true and false contacts. 


®@ Tactical platform objects which simulate the behavior of their real world 
counterparts. These are dynamic associations of component objects which may 
change in response to simulation events. A submarine inheriting the properties 
of a Patroller object may switch to an Intercepter after making a detection. 

The program flow from object to object is not sequential, but event driven. 
Once started, the simulation clock advances time relative to the simulation. Sensor 
objects begin glimpsing and motion platforms begin moving. Possible detection 
events are sent to the event manager for proper sequencing and resolution. Detections 
are communicated up the designated chain of command to data-fusion centers where 
an instance of the ACT processes them and integrates them into the tactical picture. 
This may cause a command object to direct (or redirect) assets in response. This 
process continues until the allotted time for the scenario expires. The MOE Statistics 
for that run are added to the MOE file and the next iteration of the scenario 


commences. When the desired number of iterations are complete the compiled 


summary of MOE statistics is analyzed to interpret the results of the simulation. 


2. ASW System Architecture Evaluation. 

In order to evaluate a desired architecture it must be conceptualized in great 
detail. The classic computer maxim “garbage in, garbage out" is especially true of 
ASSET where a single bad parameter value can render the results meaningless. 
While the user-friendly interface facilitates the mechanical process of inputing data, 
the abstract nature of the required parameters make it necessary to assume values of 


dubious validity at times. General areas which require quantified data include: 


Communication connectivities which represent how detection reports are passed 
from object to object and what delays are involved at each node. 


Command organization including geographic regions of responsibility. 
Environmental data for acoustic propagation loss as a function of range and 
frequency as well as ambient acoustic noise in each distinct environmental 


region. 


Motion plans representing the operating areas and interconnecting tracks which 
will govern tactical platform movement. 


Umpire parameters such as the kill probabilities for a given submarine class 


against each possible target class both for the case where it detects the enemy 
first and the case where it is detected by the enemy first. 


The individual platforms require complex definition as well. A submarine 


object, for example, requires parameters for: 


breakoff speed between "fast" and "slow" acoustic behavior. 

self-noise at slow and fast speed. 

directivity index and recognition differential (for Sonar Equation computations). 
signature frequency emitted and intensity at slow and fast speed. 


movement speeds when patrolling and intercepting and detailed motion plan to 
be followed. 


weapons loadout and level at which the sub will abort its mission to rearm. 


whether the sub will transmit the detections it makes and risk detection itself or 
not report any detections. 


the interval at which it copies the submarine broadcast for orders to investigate 
cues. [Ref. 2:p. 2-37] 


The submarine’s sensor requires a detailed set of parameters of its own which 
will be used to construct the contact report that is communicated to the ACT 
representing the data fusion center associated with the submarine’s command object: 

® position error which represents the 2 standard deviation radius of a circular 
normal distribution of contact errors about the actual target position. 
® course and speed errors which are uniform about the true value. 


®@ Target Motion Analysis (TMA) delay representing the time required to acquire 
course and speed information. 


@ false alarm rate. 
@ determination whether, on the average the sensor is reliable enough for the 
fusion center to automatically start a new track based on a single contact report. 
The complexity of the data involved in constructing an architecture for 
evaluation varies considerably, as can be seen from these examples. Once all this 
information is input, it can be easily double checked by accessing the appropriate 
windows again. When the user is sure al] data is input correctly, the simulation can 
be run with the platform graphics turned on to ensure it progresses properly. These 
graphics can then be turned off to run multiple iterations faster. 
When the desired number of iterations is completed a dialog box opens asking 
if the MOE’s are to be saved. Once saved to disk, they can be opened and examined. 


They include: 


® Attrition of submarines and MPA assets for each side as a fraction of original 
force. 


® Number of times submarines approached to within a critical range (assigned by 
the user) of enemy surface formations (which serve only as targets in ASSET 

1.0, they have no combat or detection capability). 

® Tracker statistics involving the fusion delay involved in constructing target 
tracks. 

These statistics along with weapon expenditures allow the completed simulation 
to be compared to other runs to observe the effect of a particular parameter changing, 
or how close a particular MOE comes to a goal value. The methodology of ASSET 
must be taken into account as the computer commanders follow simple resource 
allocation rules which do not necessarily reflect the tactical priorities the user desires. 
A higher priority for prosecution may well be given to a submarine returning to base 
than one ten miles from a surface group, due to the resource allocation process’ 
emphasis on detection rather than on protection. While creativity is needed in the 


design, implementation and interpretation of an ASSET scenario, significant insights 


into the conduct of ASW campaigns can be gleaned from it. 


3. Correlation, Tracking and Data Fusion. 

ASSET is designed to simulate the data fusion centers where contact reports 
from a wide variety of sensors and platforms are centrally processed to create a 
tactical picture for the region of interest. These data fusion centers are simulated by 
an instance of the OBU-ACT module with appropriate communication connectivities. 


The result of correlating and processing the streams of contact information arriving at 
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a fusion center forms the basis for the allocation and cuing of assets. Thus it is 
critical that the ACT is not spinning its wheels performing unnecessary computations. 
The implementation of OBU-ACT used in ASSET 1s organized into four main 
functional areas: 
® generating a new track from one or more contact reports (the START and 
CLUSTER modules). 


® associating a contact report to an existing track (the INPUT, START, and 
ASSOC modules). 


® simulating the intervention of a human analyst to resolve ambiguous contact-to- 
track associations (the ANALYST module). 


® updating and managing a database of the status of all contact reports and tracks 

(the Locational Data Base Manager and MTST modules). 

The data flow through these modules is shown in Figure 1. A detection is made 
by a sensor and a report is initiated. After the applicable delays have elapsed, the 
contact report travels along the communications network to the fusion center, 
experiencing additional delay at each node. 

Upon arrival at the fusion center the INPUT module enters it into the Locational 
Data Base Manager (LDBM) module and passes it to the ASSOC module. There the 
contact is compared to the existing tracks on the basis of geofeasibility. Those tracks 
with which the contact could be associated are evaluated by a statistical comparison of 
a measure of correlation (MOC) with a preset threshold. If the contact fails to trigger 
an association with any track it is passed to the START module. If the contact 


triggers a single association, or if one track MOC exceeds all others by a preset 
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Figure 2: ASSET contact data flow. 


margin, it is unambiguously associated with that track [Ref. 4:p. 3]. Ambiguous 
association is resolved by the ANALYST module, which has a given probability of 
making the correct association (.7) and an associated exponentially distributed delay 
time spent making the decision (mean of .2 hours) [Ref. 3:p. 11; Ref.5:ACT- 


ANALYST]. 
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The contacts which are unresolved are passed to START which checks the value 
of the flag single-report-to-track. If the value is true, START initiates a new track 
based on the contact, if it is false, the contact is passed to the multiple contact track 
initiation module CLUSTER. This represents a simplification of the actual OBU- 
ACT START module which uses comparison to a set of criteria and human analyst 
interaction to determine if a new track is to initiated [Ref. 6:p. 3-4]. Within the 
context of the simulation, this simplification makes sense as two classes of detecting 
sensors are represented; those that report continuous contact observations such as 
SOSUS or other fixed area sensors, which are processed at the fusion center; and 
those which perform platform-level analysis and periodically report the tracks they 
hold, such as submarines. Platforms in the first category, unless possessing 
unusually low false alarm rates, would not initiate a track based on a single report, 
while those in the second, would. 

Unassociated contacts which remain after passing through the ASSOC and 
START modules are processed by the CLUSTER module. Here subsets of all 
unassociated contacts are evaluated for geofeasibility and constant course and speed 
likelihood. If the value of the constant course and speed likelihood exceeds a 
threshold value then the subset of unassociated contacts is used to initiate a new track 
[Ref. 7:p. 1-3]. Contacts which are not associated after passing through ASSOC and 
CLUSTER carry counters which increment each time they pass through the two 
modules. When these counts exceed maximum limits (hardwired at 7 in CLUSTER 


and at 5 in ASSOC) they are dropped from the LDBM [Ref. 5: OBU-ACT]. 
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When a new track is initiated or a contact is associated with a track, the MTST 
module is called to filter the track data and produce an optimal estimate of the target’s 
location and velocity along with the error covariance matrix describing the quality of 
the estimate. The full capability tracker in MTST contains the matrix versions of the 
Kalman filter equations. These equations determine the mean target position at a 
future point in time and the covariance matrix representing the error in that 
prediction. The elements of this covariance matrix can be recast as ellipses 
representing two o or 86 percent containment regions for the target’s possible position 
and velocity. A complete description of MTST is the subject of Chapter II. 

The foregoing process of correlating contacts to tracks and filtering those tracks 
to provide an estimate of target position occurs at each instance of a data fusion 
center object. When formulating an architecture, the fusion center’s role must be 
carefully designed. Since no contact reports flow out of a fusion center, it is not 
possible to directly fuse the pictures at two or more fusion sites into a higher echelon 
center. The individual sensors originating the detection reports must send duplicate 
reports to all fusion centers which will be using the report in determining a tactical 
picture. The fusion center also keeps separate instances of the tracker to process 
surface contacts and subsurface contacts. Together with the START processing of 
single contact to track decisions, this assumes that a great deal of contact processing 
is going on below the level of the fusion center. 

Assuming that reports arriving at the fusion center are pre-processed justifies 


several assumptions that the DETECTION-REPORTER object makes about the 
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reports it generates. Most important to streamlining filter operation is the assumption 
that all submarine probability areas (SPA’s) are circular. While this represents a 
considerable inaccuracy for single contact line-of-bearing detections, if one assumes 
the contact report represents the accumulation of a significant number detections, the 
circular error assumption is more reasonable. Several types of detections are resolved 
in ASSET without explicitly modelling the geometry of the searcher and target at the 
time of detection (MPA and Fixed Area Sensors specifically). The calculation of the 
size and orientation of explicit AOU ellipses where possible would slow execution and 
add considerable complexity to the program, requiring even more specific user input 
sensor parameters. 

The assumption that this circular AOU size is independent of range is more 
difficult to defend, in light of an expected linear relationship between the standard 
deviation of the position error and range. That assumption on ASSET’s part can be 
remedied easily, for the cases where range 1s explicitly determined, by expanding the 
radius of a sensor’s AOU linearly with range based on a user input bearing error. 
This would make the process of tabulating preprocessed variance data impossible, 
however, as the sensor AOUs would no longer be constant. A good compromise 
involves defining three AOU sizes for contacts at short, medium and convergence 
zone ranges. This would allow for some representation of the range effect on AOU 
size, with acceptable complication of the variance pre-calculation procedure. 

The other tracking assumption ASSET makes is that the covariance of the 


velocity can be represented by the steady-state stationary cross-covariance of the IOU 
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process, independent of the sensor characteristics. This assumption is used to 
initialize the covariance matrix of a new track. This method is correct for a position- 
only contact which gives no information about velocity. Given position and velocity 
information however, this is only accurate if the sensor is so inaccurate that the IOU 
process velocity variance is small compared to the sensor velocity variance [Ref. 8:p. 
5-20]. 

This is not true in most cases however, with the result — that the filter re- 
initializes the velocity update each time a contact is processed, never including any 
representation of the velocity accuracy of the sensor involved. Thus, the velocity 
variance is based solely on the IOU model and not on the variance of the velocity 
reported by the sensors. This assumption makes the tracker more responsive to 
course changes, but it also artificially inflates the SPA size of accurate contact 
streams because it ignores the sensors velocity error, which may be significantly less 
than the JOU value. 

Given these assumptions and the noise characteristics of the IOU process, the 
filter implementation can be streamlined to better match the data provided. This 
process is detailed in the following chapters which will describe the IOU process and 


the formulation of the general and ASSET-specific filter equations. 
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Il. THE MANEUVERING TARGET STATISTICAL TRACKER 


Maneuvering Target Statistical Tracker (MTST) is the name given to a class of 
Kalman filter-based estimation and prediction algorithms utilizing the Integrated 
Ornstein-Uhlenbeck process to model the statistics of target motion. To successfully 
track a maneuvering target, the details of the target path should be statistically 
predictable. How accurately an automated tracker like MTST performs its task is 
directly related to how well the tracker’s target motion model describes the targets 
actual random Deca 

The random tour model [Ref. 9] can be used to describe the statistical properties 
of a target moving at constant speed which makes random course changes at times 
separated by intervals chosen from an exponential distribution. This provides a 
reasonable estimate of the type of motion expected of a patrolling submarine target. 
Unfortunately, the target distribution generated by a random tour is not normal and 
cannot be used directly with a Kalman filter [Ref. 8:p. 2-10]. A process is required 
which produces normal distributions which best approximate the mean and variance of 
the random tour. This is the IOU process. The mathematical development which 
follows is summarized from the more complete treatment in Reference 8, the results 


of which are in agreement with References 5, 10 and 13. 
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A. THE INTEGRATED ORNSTEIN-UHLENBECK PROCESS 


The Integrated Ornstein-Uhlenbeck Process (IOU) is based on a first order 
stochastic linear differential equation which approximates the higher order non-linear 
deferential equation which defines the motion of a randomly maneuvering target. It is 


a member of a class of equations known as Langevin equations and has the form: 
Ful = -Bu(t)dt + o dw(t). (2.1) 


The first term represents a deceleration caused by a resistive force proportional to the 
velocity u(t). The random forces acting on the particle are represented by the second 
term, where w(t) is a Gaussian white noise process. 

The stochastic nature of this equation makes it possible to find only the statistics 
of the distribution of the solutions, rather than the specific solution itself. A Gaussian 
w(t) produces a Gaussian velocity process u(t) which, like any normal process, is 
defined completely by its first and second order moments. The result is the velocity 
distribution of a particle which is undergoing random motion similar to Brownian 
motion, experiencing random instantaneous accelerations, but whose velocity is 
damped by a spring-like effect which constantly accelerates the particle in the 
direction opposite its velocity at a rate proportional to that velocity. The position 
variance of such a particle is unbounded over time, but the velocity variance is 
limited by the damping coefficient 8. The result is a velocity distribution that is 


normal with a limiting variance given by: 
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lim V{u()} = limE {{u(@)-E {u@}P} = (2.2) 


2B 
and a mean for any given ¢ given by: 
E{u()} = e“E {u(0)}. (2.3) 
The mean and variance expressions for the position of an object whose motion 
is governed by an IOU process completely specify its positional distribution over 
time. In order to show that the IOU process approximates the same motion as a 


random tour, the radial position distributions for the two stochastic processes will be 


shown to have the same variance. 
B. EQUIVALENCE OF THE RANDOM TOUR AND IOU PROCESSES 
The variance of the radial distance from (x(0),y(0)) is: 


E{R(p} = E{x()-xOP+[y(-yOF}. ea) 


For the Random Tour, the following holds: 
E{R%0} = 22a t-1+e-*4, (2.5) 
a 


where V is the speed of the randomly touring particle and a is the mean number of 


course changes per unit time. The corresponding result for the IOU process is: 


19 


E{R(n)} = 18 t-1+e FF}, (2.6) 


where o is the scale factor for the random acceleration and £ is the damping constant 
on the velocity, as in Equation (2.1). 
These two equations can be made identical given the following relationships 


hold: 


= V (2.7) 


| 4 


Thus, a Random Tour with parameters a and V can be approximated by an IOU 
process with parameters 8 = a and o = Va'”. It is worth noting that since a normal 
distribution is completely specified by its mean and variance, the IOU process 
represents the best normal approximation to the Random Tour [Ref. 10). 

In ASSET the IOU parameters, 8 and o are hardwired to represent a target 
conducting a random tour at ten knots with an average time between course changes 
of four hours. These parameters should match the actual motion of the platforms 
modelled in a given scenario. Since the user can choose whatever target motion 
parameters he wants for each individual platform, the parameters embedded in 
ASSET may disagree considerably with the actual target motion. This disagreement 
causes the IOU prediction to consistently lead or lag the target’s mean position, 
depending on whether the model’s speed and course change rate correspond to a 
velocity faster or slower than the target’s speed, respectively. This causes the 


position of the center of the resulting SPA to be offset from the true target position. 
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Even when the tracker parameters exactly match the target motion, a considerable lag 
occurs if the time between measurements is not significantly greater than the time 


between maneuvers. 


C. OTHER MOTION MODEL OPTIONS. 


The IOU model used by MTST damps out the mean speed of the target 
exponentially over time to account for target maneuvers. No matter how small the 
time interval, the process reduces the mean speed by the proper amount. If the 
interval between detections is consistently shorter than the maneuver interval, the 
damped velocity is used to predict the position of the constant velocity target and the 
tracker develops a lag. One remedy for these inconsistencies would involve 
increasing the complexity of the tracker by using adaptive methods to more closely 
match the model parameters to a given target’s track. An adaptive filter recognizes 
changes in the targets motion and compensates for it by changing the parameters of 
the motion model. 

Desiring to reduce the filter’s complexity however, instead of explicitly 
modelling an adaptive tracker, the use of such a tracker can be approximated using 
the IOU process as a basis for the size of the SPA, but assuming the model has some 
adaptive properties to position it more accurately. This can be accomplished using 
the existing IOU model to estimate the SPA variance and by assuming the tracker’s 
position prediction’s are distributed about the target’s true position in a circular 


normal fashion. This is an optimistic assumption, but captures the fact that real world 
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trackers, using more complete contact data than ASSET and an adaptive motion 
model, can achieve more accurate position and velocity estimates. Examinations of 
such trackers can be found in References 14, 15 and 16. Reference 14 in particular 
demonstrates a passive tracker that delivers outstanding accuracy while the contact is 
on a constant course and speed, but has difficulty maintaining track during difficult 
target maneuvers. 

The benefits of adaptive tracking are not applied to the SPA size in this case, as 
the ASSET tracker is not an adaptive one, but are assumed in simplifying how its 
position is determined. This technique introduces a moderate degree of error in SPA 
size by not precisely simulating the mechanics of distributing the position of the 
center of the SPA. It also does not take into account the filter’s velocity components. 
Reference 15 provides methods for countering these effects through noise adaptation 
and correlated maneuver gating while Reference 16 uses bias-sensitive maneuver 
detection and Kalman gain adaptation. Another option used to improve the IOU 
model is the use of a dual velocity system which combines a short term IOU process 
combined with a long term one. This type of model can reproduce a variety of 
motions depending on the weighting factors given to the short and long term 
components. The resultant Kalman filter has six states and a 6 x 6 covariance matrix. 
This is the implementation of the IOU process currently used in the TOMAHAWK 
fire control system for surface ship motion modelling [Ref. 8]. These are but a few 
of the schemes more complex trackers can use to improve motion modelling. The 


single velocity IOU process, while not adaptive, comes close to satisfying the 


Ze 


necessary conditions and since it is the basis of the ASSET tracker will be used for 


comparison. 
D. THE KALMAN FILTER 


In order to estimate target position, a method of combining the contact 
position and its covariance with a prediction of that position and covariance based on 
the IOU process is necessary. The Kalman Filter provides such an estimate, provided 
the process model is linear and time invariant, which the IOU process is. The 
Kalman filter is derived from the least-squares estimation theory of Gauss and 
Legendre [Ref. 11:p. 7]. This technique sought the most probable value of an 
unknown quantity based on a series of measurements containing unknown errors. 


More formally, finding the most probable estimate of x based on the matrix equation: 
z, = HAX+Y, (2.8) 


where vector Z, represents a measurement taken at a discrete time kX, based on an 
observation matrix H revealing information about one or more state parameters of the 
true state vector x but corrupted by a noise vector v. This technique was generalized 
by Kalman to apply to linear filter theory, producing a filter which produced a best 
linear estimate of x [Ref. 12]. 

The resulting filter is based on two mechanisms operating sequentially, one 
predicting the state and covariance for a future time and the other combining this 


prediction with a noisy measurement of the target’s state taken at that time. This 
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should, given correct motion modelling and noise parameters, produce a state and 
covariance which more accurately reflects the true target state than either the 
prediction or the measurement alone. 

These mechanisms are represented by two systems of matrix equations. The 
variables involved are defined in Table 1. The prediction equations for the state x 


and covariance »& are: 


state prediction: %,,, = 6X,+H,, (29) 
covariance prediction: ,,, = ¢2,¢7+Q,, (2.10) 


where ¢ represents the transition matrix from the dynamic system model, Q 
represents the variance inherent in the dynamic model and yz,, represents the mean of 
the noise term w from the dynamic model, which is assumed to be zero in ASSET. 
The equations for combining the predicted and measured states and variances require 
the computation of a Kalman gain matrix which provides the weighting factors for the 


combination. These equations which update the filter’s estimates are: 


A A -1 
Kalman Gain: K,,, = 4,,,4 lS, HT +R, (2.11) 


(2.12) 
State update: x,,, = %,.,+K,., [Zeet “i, - HZ,,,], 


covariance update: ¥,,, = [I - K,.,]¥,.1, | (2.13) 


where KR represents the variance of the noise inherent in the measurements, Z is the 


measurement itself, 4, is the mean of the noise inherent in the measurement (also 
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Table I: Kalman Filter Definition of Terms 


Kalman Term: General Form ASSET Form 
State vector at time k: tei ke! mean 4 x 1 


State transition matrix, 
time k to time k+/: @ nxn phi 4x4 


State process noise mean: wy nx 1 assumed zero 


State error covariance, 
at time k: yUnxn var 4x4 


State process error 


covariance at time k: Onxn f4x4 
Observation matrix: Hmxn I 4x4 
Observation at time k: Zmx 1 state 4x 1 
Observation error mean: py, mx 1 assumed zero 


Observation error covariance 
at time: Romxn covariance 4 x 4 





assumed to be zero in ASSET) and J is the identity matrix of proper size. The state 
and covariance updates, assuming H (as ASSET does) is the appropriate identity 


matrix can be written as: 


state update: x,,, = 7 el Pee ae + Rua "Zpa): (2.14) 


A - -1 
covariance update: %,,, = (5. +R, 7) (2.15) 


where the inverse of R,,, and L,,, are the weighting matrices. Thus the inverse of 


the noise error variances are the weighting factors by which the prediction and 
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measurement are linearly combined. This is the form of the Kalman equations which 
ASSET actually uses. The advantage of this form is the ready understanding of the 
underlying mechanism of weighted linear combination. The disadvantage is the 
computational necessity of performing three matrix inversions to calculate the 
variance propagation. Regardless of the form of the equations, the ¢, QandR 


matrices must be developed to complete the computational form of the filter. 


1. Development of the ¢ Matrix. 

The matrix ¢ represents the transition matrix which governs the dynamics of the 
state between two discrete times. The following development summarizes the 
derivation of Reference 11, Section 5.2. The differential equations defining the IOU 
process and velocity can be combined to form a single matrix differential equation 


which takes the following form in the single coordinate direction x: 


0 ] 
a _ | , "| 58 2.16) 
BLU) 0 8B u(t) o 


more compactly: X(t) = F(t) X(t) + Gwi(t) (2.17) 
where X(t)represents the system state vector and the dot over a symbol indicates 


differentiation with respect to time. Given that X(t) has the value X(¢,) at time ¢, then 


X(t) at any future time f>12, is given by: 
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X(t) = (t,t) X(t) + | $(t,7) Gw(r) dr (2.18) 


Th, 
where $(t,t,) represents the transition matrix between the states at time ¢, and time ¢. 
Provided that F(t) is constant in time and assuming that f, is zero, the matrix ¢ can be 
shown [REF. 8:p. 5-3,5-5] to have the form: 


1 1ta-e-# 
B (2.19) 


This represents the transformation projecting a state into the future based on the IOU 
model. The velocity mean decays to a new value which is simply the old value 
multiplied by e?'. The position is moved a distance equal to the time 1/8 multiplied 


by the change in velocity: 


= 


3 Me MOeo) (2.20) 


Xx(0),., = x(t),+ 
where: 


u(t),,, = u(,e*. (2.21) 


The form this takes when expanded to cover both coordinate directions is: 


Vii 


B 
1 c 
o(t) = 0 ] 0 em Br) : (2.22) 
0 0 exe 0 
0 0 0 is 


where ¢ is time elapsed from the last measurement to the observation time of the 
current measurement. It is important to notice that the transition matrix depends only 
on the IOU parameter 6 and not on o. Thus the IOU velocity damping coefficient 
and the elapsed time determine the dynamics of the mean state and variance 
transition. Appendix A contains graphs of the values ¢(1 2), which will be called ¢2, 
and $(2),which will be called $3, take over various time intervals. The noise 


magnitude coefficient o is involved in the noise terms R and Q. 


2. Development of the @ Matrix. 

The Q matrix represents the noise introduced by the white Gaussian process 
w(t). As discussed above, as the elements of this matrix get large the filter 
increasingly relies on the measurement data and ignores its predictions. As the 
elements get small the filter increasingly ignores the new contact data and the errors 
in its predictions compound until it diverges, completely losing track of where its 


target is. The definition of this matrix, the process noise error covariance, is: 
ii VE F 3} 
0 = El@-E@} {z-E@)'| te 


which for the random process operating in (2.26) becomes: 
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OVE [serene | susp Owe | (2.24) 


T=l, S*i, 


Since the mean of the white Gaussian processes w(t) and w(s) are both zero and the 


variables of integration 7 and s are equal, this reduces to the form: 


Q = | $(t,7r) GGT o"\t,r) dr. (2.25) 


Substituting the appropriate matrices G and ¢ from (2.25) and (2.29) respectively and 


performing the indicated matrix multiplication results in: 


Or d -@ Tih Lar -e “A @ -B(t-7) 
B f 


t p? 
Q = J dr 


on l-e “BEM @ ~B(i-7) ore ~2BU-n 


(2.26) 


for the single coordinate direction x. Assuming that ¢,=0 and recognizing that q(1,2) 
is equal to q(2,1) results in three expressions to integrate to get the final value of the 


Q matrix in one direction: 


qll-=—s q12 
Oe where ql2 = q21 (2.27) 
q2l = =q22 
git = 2) 21-1 -4e-# +e (2.28) 
2p? B 
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o 


q12 = q21 = oa! - Qe -B + @ 261) (2.29) 
= o — e ~2ht 2.30 
q22 xg e ~2h1) (2.30) 


When projected into two coordinate directions the full representation of the Q matrix 


1S: 


Q = (2.31) 


Each element in this matrix is directly proportional to o? and decreases with 
increasing 8. Time has an inverse exponential relationship causing the noise to 
increase the longer the interval between measurements. At ¢ equal to zero the noise 


terms are equal to zero and as ¢ goes to infinity : 


o o? 
ql]1-~, ee q22-— (252) 


As mentioned previously, the variance in position is unbounded, but the cross- 
covariance and velocity covariance are bounded. These values become important in 
developing the measurement noise covariance matrix R. Appendix A contains graphs 


of these values over various time intervals. 
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3. Development of the R Matrix and Filter Initialization. 

The matrix R represents the covariance of the measurement noise and 
accompanies each state measurement reflecting its precision. This value may vary 
with time in the ASSET tracker as contact reports flow in from a variety of sensors. 
The inverse of this matrix is used as the weighting function for the contact report, the 
larger the elements of R the less effect the contact has on the filtered update. The R 
matrix also plays a role in initializing the filter. 

When a target is first detected, there is no f-J state for the filter to work from in 
making its prediction. In general, a time ¢, R matrix is preset to initiate the filter. 
This matrix would consist of the 2 x 2 positional error matrix which describes an 


estimate of the elliptical AOU about the contact position: 


1 “ cos’6) + (B* sin’?@) (B* - A’) sin@ cosé | (2.33) 
4 | (B?-A?)sin@cos@ (A? sin26) + (B? cos’6) 

For a position-only contact, the two lower right diagonal entries would be set equal to 
the stationary variance of the IOU process (2.2). This represents the maximum 
variance the velocity could possibly have based on the motion model and since no 
velocity information is coming in, the filter is relying solely on the motion model, 


making this a value a reasonable choice: 
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EO, 203 le 
—_ 0 0 
(2.34) 

o 

is; 

0 0 = 
ow 
2 
alto - 


If the contact is a position-and-velocity contact, the procedure is more complex. In 
like fashion to the initialization of the upper left 2 x 2 portion with the variance in the 
sensor (its elliptical AOU), so must the bottom nght 2 x 2 be initialized by the 
velocity variance of the sensor. Since the IOU process is predicting the velocity 
distribution also, this variance of the measurement must be combined with the 
variance of the IOU process. The proper combination can be shown [Ref. 8:App. F] 
to be the same as the alternate form of the Kalman covariance update equation (2.14). 
This method is not used by the tracker in ASSET. The ASSET tracker simply 
uses the initial contact’s IOU velocity variance as the basis for a new track, without 
updating it. This original contact report is used as the f-J state and covariance when 
the next contact associated with the track comes in. The contact reports themselves 
also differ from the construction described above. The positional variance of the 
sensors are all circles in ASSET so the 2 x 2 matrix which occupies the upper left of 


the R matrix is actually the 2 x 2 in the upper left of: 
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a 0, 60 
4 
9 
0+ 0 0 
(2.35) 
o” 
Vy edie 
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o 
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The velocity variance in the lower right hand 2 x 2 is correctly set for a position-only 
contact, but this is also used for a position-and-velocity contact. Thus the R matrix 
which ASSET uses never incorporates the error in velocity inherent in the sensor. 
Using these values essentially reinitializes the filter each time a measurement is made. 
This is desirable for position-only contacts as it prevents divergence of the velocity 
components of the state, however for position and velocity contacts, this can cause 
poor velocity accuracy as the tracker improperly weights the velocity components of 
the contact report. The velocity information of different sensors are given equal 
weight regardless of the relative size of the errors in the reports they make. The 
actual velocity variance of the sensor should be computed and used to set the value of 


the filtered velocity using the form of (2.15). 
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Tt. IMPLEMENTATION OF THE TRACKER IN ASSET 


Having developed the components of the Kalman filter, the construction of the 
tracker implemented in ASSET can be discussed. The following is based on 
Reference 5, the ASSET source code for version 1.0. The ACT-MTST module 
consists of objects which perform the filter calculations described above, as well as 
objects which maintain the spatial relationships between the positions involved on a 
spherical earth. 

There are several constants which are defined for use in all the module’s 
objects. The first is dtor which is used for degree to radian conversion. Where 
trigonometric functions are indicated below, the actual code uses dtor to convert 
angles from degrees to radians, but in the interest of brevity, this will be excluded. 
Second is a, which is equivalent to 8 in the IOU process, and is set to .25. Third is 
o which is the noise scale factor from the IOU process as described above and set to 
V50. This is equivalent to the square root of a velocity multiplied by the Random 
Tour velocity. Thus the tracker predicts the target position by assuming a random 
tour 1s being conducted at a speed of ten knots with a mean interval between course 
changes of four hours. This cannot be changed by the user. 

Inputs to ACT-MTST are read from the data in ACT-LDBM, the contact 


database manager. Ifa new track is to be initiated, the next MTST object, start-new- 
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track sets the track head and track tail state and covariance matrices equal to the 
contact mean and covariance using compute-mean-covariance. The use of two sets 
of data, a head and a tail, to allow incorporation of out-of-sequence contact reports 
will be discussed below. The object compute-mean-covariance is also called by 
ACT-ASSOC as a basis for the spatial measure of correlation used to determine 
contact to track correlation. 

Compute-mean-covariance creates two arrays, a four by four called var, which 


is the R matrix, and a four by one called mean, which is z,: 


0 
0 (3.1) 


mean = | spd cos(90-cse) | ’ 
spd sin(90-cse) 


and var, which has the same form as (2.34). The vector mean is Centered at the 
latitude and longitude indicated in the contact data list and x and y velocity computed 
from the contact course and speed. This variance matrix however, as discussed 
above, is calculated as if the contact report consists of an ellipse with minor axis B, 
major axis A and orientation angle @. All the spatial data fields that are generated by 
the Detection Reporter and Generic Sensor contain a single positional uncertainty and 
an orientation angle of zero. Modifying (3.1) to take advantage of this fact results in 
(2.35) where r is equal to the standard deviation of the sensor’s target position 
estimates. This alternative form returns matrix element values identical to those 


computed by (3.1), given the actual form of the contact reports generated. 
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If a contact is associated with an existing track rather than used to start a new 
one, the update-track object is called by ACT-LDBM to perform the Kalman filter 
operations. The remaining objects in ACT-MTST are all components of the update- 
track object with the exception of the object cov-to-ellipse, which is called by various 
display and graphics objects to determine the orientation, semi-major and major axes 
of the ellipse representing a covariance matrix. Since, as shown above, all the 
ellipses are actually circles, this object is another anachronism of the mismatch 
between the capabilities of the tracker and the contact reports actually input to it. It 
can be simplified in a manner similar to compute-mean covariance above to return 
the radius of a circle instead of a major axis, minor axis and orientation of a ellipse. 

Focusing once again on the update-track object, it first defines the spatial data 
parameters of the contact and the track the contact is used to update. ASSET uses 
four data structures to represent the data associated with contacts and tracks: 

® obu-contact: consisting of contact id, receipt time, track association, sensor, 
categorization (track association flags “pinnedp" and “lockedp", and number of 
passes through ACT-CLUSTER and ACT-ASSOC), spatial data (the obu-data- 
field described below), altitude (surface or subsurface) and HFDFp (an optional 
parameter associated with HFDF contacts). 

® obu-data-field: type (position-only or position-and-velocity), observation time, 
latitude, longitude, AOU orientation, AOU major axis, AOU minor axis, 
course, course uncertainty, speed and speed uncertainty. 

® obu-track: track id, number of contacts associated to the track (maximum of 
five), state and covariance of the track head (most recent contact) and track tail 


(oldest online contact), head and tail spatial data (each an obu-data-field), head 
and tail contact id, and altitude (surface or subsurface). 
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® obu-state-field: reference latitude, reference longitude, state array (offset from 
reference position in miles N/S, offset from reference position in miles E/W, 
N/S velocity, E/W velocity) and covariance array (four by four array per (3.1)). 
The three AOU variables can be eliminated from obu-data-field and replaced by 
position-uncertainty, the radius of the actual AOU circle reported. This will not 
affect the tracker’s calculation and reduces the memory taken up by contact reports. 
An obu-track consists of up to five contact reports, the most recent of which is 
the track head and the oldest is the track tail. The variables required for performing 
the filter calculations are read from the data structures described above, update-track 
then checks the time since the last update. If the new contact observation time is later 
than the track head update time, the contact 1s used to update the track head using a 
set sequence of object calls. This sequence contains the Kalman filter procedures as 
well as two routines which maintain the latitude and longitude relationships on a 
spherical earth and proceeds as follows: 
® ioumotion: Performs the IOU prediction step and updates the N/S and E/W 
offset distances in the track head state to the time of observation of the contact. 
These offsets effectively perform the prediction in a plane tangent to the earth’s 
surface at the reference latitude and longitude point of the track head. 
® change-tangent: Finds the relationship of the iou predicted position and the 
contact point in a plane tangent to the earth at the contact’s reference point vice 
the track head’s reference point used to compute the prediction in ioumotion. 
These matrices are stored as state-cov-old. Since the covariance is circular, 
there is no need for the sequence which rotates it. 
® make-state: Takes the contact data and produces z, which is called simply 


mean and its covariance R which is called var. These two matrices are stored 
in state-cov-ctc. 
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filter: Takes state-cov-old and state-cov-ctc and uses Equations 2.14 and 2.15 
to compute x,,, and 4,,,. The result is the Kalman filter update performed in a 
plane tangent to the earth at the contact’s reference point. 


center-tangent: The latitude and longitude of the position produced in filter is 
projected back onto the earth from the tangent plane. 


update-track-head: The new updated state and covariance are made the new 
track head. The contact used to perform the update provides the track-head 
contact data. 


update-track-tail: Sorts the contacts making up the track in time order, oldest 
to most recent (1 to n). If there is only a single contact, it is designated the 
track tail. If there are five or fewer contacts and the oldest contact is already 
the tail it remains so; otherwise the oldest contact is made the new tail. If there 
are more than five contacts, it takes number contacts, subtracts five and the 
contact whose number equals the result is used to update the track tail using the 
procedure detailed above. The contact used to update the tail provides the track 
tail contact data. All contacts older than the tail are then pruned from the 
track. 


Considerable computation time can be saved by simulating a portion of the above 


process, rather than actually carrying out all the steps. A revised update procedure 


based on using the updated variance (2,,,) as an approximation to the distribution of 


updated states (x,,,) 18 proposed. Rather than produce the contact position latitude 


and longitude randomly using the sensor variance as a distribution and processing it 


using the sequence above, the latitude and longitude are chosen by computing the 


updated variance, which can be done prior to scenario execution. When it is needed, 


the distribution is centered on the target’s true position, and used to generate the 


updated states directly by random draw. This new procedure analyzed in detail in 


Chapter IV. 
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Time delays incurred at communications nodes enroute to the fusion center may 
cause the observation time to fall between the track head update time and the track 
tail update time. In that case a modified procedure is followed . 

® If the time of observation for the contact is later than the tail update time or if 
the track does not yet consist of five contacts, then the contact is appended to 
the track and the contacts are sorted by time of observation. If the contact is 
older than the tail and the track already consists of five contacts, it is ignored 
and the procedure stops. 

@ If the observation time is later than or equal to the tail update time, the tail is 
used as the oldest contact, if the new contact is older, then compute-mean-var 
is used to make the new contact the oldest. 

@ The oldest contact is then updated using the procedure described above using the 
next contact in the track, in time sequence, as the "contact". The result of this 
becomes the “old contact" and the procedure is repeated, updating each contact 
in the track in sequence until the head is updated. The resulting "forward 
filtered" track closely approximates the track head which would have resulted if 
the time-late contact had been received in sequence. 

The calculations performed in ioumotion and filter, all assume that elliptical 
AOU’s are being produced by the Detection Reporter and General Sensor objects. 
The object make-state can be simplified in the manner as compute-mean-cov as it 


performs the same operations. The reduced forms of the equations in those objects 


are computed below. 


A. REDUCED COMPLEXITY IOU PREDICTION 


The matrix equations for calculating the IOU predicted state and covariance are 
given by Equations 2.9 and 2.10. Using the symbolic processing capability of the 


MathCad 3.0 (copyright 1991 MathSoft, Cambridge MA) computer software package, 
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the matrix operations were evaluated. The large number of zero terms and high 
symmetry of the ¢ and Q matrices resulted in simple expressions for the elements of 
the state and variance predictions. Portions of the calc-mot-mat object which 
compute f(0,0), f(0,2), f(2,2) for positive time as well as ¢(0,2) and $¢(2,2) can be 
used to produce the terms of the following equations used to compute new-var: 


E/W PoOSs,,, Vest y, +phi(O,2) Vy, 
N/S pos,,, Lia x,+phi(0,2) Vx, (3.2) 


F=newmean=! erw vel |~|Wy,.| | Vy,phi(2,2) |? 
N/S vel,,, Vy, 01 Vx, Dhi(2,2) 


where: 
4(0,2) = qare™ and (2,2) = e-* (3.3) 
and: 
Va,., 0 Vb. 0 
y. =newvar= in nadie (3.4) 
vil Vb,,, O Ve,,, 0 
0 Vb,.1 0 Vey. 
where: 
Va,., Va, +2 (phi(0,2)) Vb,+phi(0,2)? Vc, +f(0,0) 
Vb,.,| = (Vb, +phi(0,2) Vc,) phi(2,2)+f(0,2) (er) 
Ve,.1 phi(2,2) Ve,+f(2,2) 
and with: 
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|p Or ah oe ee 
iO, 0) py pie” Bme |. 
o 


= 2 | (1 _¢-ey4(1e- 3.6, 

f0,2) = Ge "Y+(se ")) (3.6) 
- © 4 2-28 

ft2,2) 73 Cue): 


The f terms given here are the ASSET versions of (2.36), (2.37) and (2.38). 
B. REDUCED COMPLEXITY KALMAN UPDATE 


In a manner analogous to that used in the preceding treatment of the IOU 


update, the matrix equations for the Kalman update were evaluated. The result is: 


Ka,,, O Kobl,,, 0 

0 Ka,,, 0 Kbl,,, (3. 7) 
Koo Ome Ke. = QT 

O Kb2,,, 0 Ke¢,,, 


Kalman Gain = 


where: 
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Ka,,, =s Va. He FM ’ 


D me B 
Kbl,., = -< Vosag P|. 
Kb2,,, -- VS | (3.8) 
Key. 3(¥ Ss Want 529M , 
D =(Va,,,+ a? 2\(Ve,..+ : 7)" Vb... 
Using these gain values the corresponding updated variance terms are 
Va,,, = Ka,,, — | 2 
VbI,., = Kbl,, £ 
LK (3.9) 
Vb2,,, = Kb2,,, ar 


x 
Vey zi K¢,., 2B 


The updated state is written in terms of residuals which are the differences between 


the contact states and predicted states: 
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Ry zy -Xy 


X-KX 
Residuals = = : 
RVy ZVy-XVy 
RVx| | zVx-£Vx 
§+Ka Ry+Kb1 RVy 
x +Ka Rx+Kb1 RVx 
Vy| | Wy+Kb2 Ry+Kc RVy 
Vx | Vx+Kb2 Rx+Ke RVx | 


(3.10) 


(Sa) 


The reduced form equations for the state are included here, but are unnecessary if the 


estimation of the state is to be used. The accuracy of the estimation technique is 


analyzed in the following chapter. 
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IV. MODIFICATIONS TO THE ASSET TRACKER 


The tracker currently used in ASSET makes two major assumptions relating to 


detection modelling: 


@ Sensor AOU’s are circles vice ellipses, 


@ The velocity variance inherent in the sensor is ignored when making the filtered 
variance calculation for position and velocity contacts. 


Figure 2 shows the geometry of a typical filter operation: 


Contact Position 
Randomly Orawn From 


Dietribution Centered on xX. 
Target’e True Poeition 


New Filtered Poeition/ 


and SPA Predicted Position ~ 


and AOU en 
Contact Position 


ae 
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Figure 2: Geometry of Filter Update Operation. 


The equations which exactly determine the position of the center of the SPA are given 
in (3.12). Given that the simulation knows the true target location and making 


several assumptions detailed below, it is proposed to bypass the Kalman state update 
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step and estimate one possible position of the SPA based on the relationship between 
the SPA center and the true target position. If the true position is distributed 
somewhere in the SPA, then from the point of view of the target, the SPA center is 
likewise distributed about its position. This allows an estimate of a SPA center to be 
generated by drawing a random variable from a distribution with standard deviation 
equal to the square root of the filtered variance in position with a mean centered on 
the target’s true position. This uses the SPA "in reverse" to find a possible position 
for its own center based on the targets true position. Provided the necessary 
assumptions are met, this will provide a reasonable approximation to the filtered track 
over a Series of observations. It also allows the variance values to be precomputed 
and drawn from a table vice calculated when needed. These three assumptions are: 

@ The mean predicted position and velocity must be very close to the target’s true 
position and velocity so the AOU of the predicted position is concentric with the 
contact AOU, or nearly so. This assumption relates to the motion model in the 
tracker simulated. 

® It is sufficient to model the behavior of detections in the mean, rather than 
explicitly computing them individually. This assumption relates to the 
interaction between the correlator and the tracker. 

® The distribution of filtered states has its mean at the target’s true position and a 
variance equal to the steady-state value of Va (3.10). 

The first assumption sets a requirement for the accuracy of the motion model. 


If the predicted position is consistently offset a considerable distance from the target’s 


true position and the Kalman gain is small, then the center of the filtered position 
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distribution will be drawn off the target’s true position. This can be shown by 


rearranging the terms in Equation (3.12) to yield for a single dimension: 


y = ~(1-Ka) + Ka Zy + Kb1 (ZVy-Wy) (4.1) 
The filtered position y approaches the contact position Zy as Ka goes to one and 
moves off to the predicted position f as Ka goes to zero. This case must be avoided 
to prevent significant inconsistencies between the distribution of actual filtered 
positions and the estimates generated. Working to keep this relationship is the effect 
shown in (2.5) which acts to push the position toward the observed position as the 
time interval increases for a given observation AOU size and the IOU prediction 
variance gets large compared to it. 

This effect is not easily identifiable in the Kalman gain equations, but increasing 
the time interval increases the size of the predicted AOU and increases Ka for a given 
observation AOU size. This counteracts the tendency of the predicted position to 
recede to the previous filtered position as the time interval increases and keeps the 
required relationship intact for cases where the time interval between observations is 
large. The maximum error will occur when the two AOU’s are nearly the same size 
and the value of Ka is approximately .5-.65. This assumption is also involved in the 
computation of the filtered velocity, where it is assumed that the difference between 
the predicted and observed position multiplied by gain Kb2 is small. 

The second assumption requires acceptance of a lack of consistency between the 


correlator and the tracker models. The ACT-ASSOC module uses the Kalman filter 
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for calculation of the spatial measure of correlation it uses to associate contacts to 
tracks. The effect of using equally representative, but different, implied contact 
locations by the correlator and tracker to do contact association and track updating has 
not been thoroughly examined. The fact that the filtered position produced by this 
method does not relate to the contact position used for track association does not 
appear to be a serious problem, given the Monte Carlo nature of the simulation. The 
draw used by the correlator can be looked at as determining whether the proper track 
assignment is made, a false assignment made or no assignment made. If the contact 
is drawn from the proper distribution, the probability that any single draw results in 
any of these occurrences is the same whether or not the tracker goes on to use the 
contact to update the track chosen. 

The third assumption is again linked to the accuracy of the motion model. The 
result of making this assumption is shown in Figure 3. The first part is essentially 
the effect of assumption one being true. The second part is required because the 
variance in distribution of the actual filtered position is affected by the speed of the 
actual target. For target speeds less than the modelled speed, the actual distribution is 
tighter than the estimated value. The opposite is true if the target speed is faster than 
the model. The implication of this assumption is that the target is exactly performing 
the motion modelled and the tracker does not respond as the actual Kalman filter does 
to deviations from the expected motion. Since the range of velocities the submarine 
targets may have is small, this deviation is typically less than fifteen percent of the 


SPA size. 
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Figure 3: Modified Filter Operation Geometry. 


A. PROPOSED MODIFICATIONS 


Hitiitementiiln of these changes in the code is a part of a separate effort to port 
ASSET to the Sun UNIX workstation. A listing of the applicable code is contained in 
Appendix B. The update-track object described in Chapter III ia the primary object 
to be modified to implement these changes. Its function can be split between those 
used to precalculate the variance values and those used to estimate the velocity if 
precalculation is desired. This involves computing the three steady state variance 


values for each sensor AOU size, for representative values of the contact interarrival 
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time. Other objects to be changed are mentioned below and included in the Appendix 


B as well. A summary of these modifications include: 


® ioumotion: The three element values of the f matrix (f(0,0), (0,2) and f(2,2)) 
and the two ¢ matrix values (¢(2,2) and ¢(0,2)) are computed as before, but not 
set in matrix form. The values are used to compute the three element values of 
the predicted variance using Equation (3.6). The three computed values will be 
used to compute the steady-state filtered variance values which are stored in a 
table for future reference. A table of the number of iterations needed to reach 
Steady state versus interarrival time is detailed in Appendix C. 


® change-tangent: This object is not required. 

® make-state: The random numbers used to determine the contact’s latitude and 
longitude (from the contactPositionDraw function called by the Detection 
Reporter) can be retained, but as discussed above this seems unnecessary. The 
two values of the R matrix are all that need be calculated by this object. 

@ filter: The values computed in ioumotion are used to calculate the three 
elements of the filtered variance using Equation 3.10 for representative values 
of the contact report interarrival time. A graph of the steady state gain Ka 
versus interarrival time is contained in Appendix D. The Kalman gains are 
computed in order to find those values and will be used below to set the range 
of operation of the filter. 

This set of routines comprise the pre-processing module and can be separated 
from the rest of MTST so it computes the necessary data at the start of the 
simulation. This requires computing the three variance values for 33 interarrival 
times, giving approximately .03 between succesive table values from .03 to .97. 
This requires 100 values be tabulated for each AOU defined. Retrieving values from 
a data array of perhaps 1500 elements appears to be significantly faster than 


calculating the values individually each time the filter is called. 
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The following routines perform the necessary computations at the time the filter 


operation is performed: 


© center-tangent: This object is modified to call the function 
contactPositionDraw and, using the filtered variance drawn from the table 
computed above, as positional variance choose a center for the SPA. The 
distances involved to not require the use of the spherical earth routine in 
contactPositionDraw. The new latitude and longitude can be calculated using 
1° of latitude=60.1077nm and 1° of longitude=60.1077*cos(latitude)nm. The 
error in this estimate is less than one percent out to ranges of 600nm at non- 
polar latitudes. The filtered velocity can be estimated by the following 
equation: 


Wie = Vy,e™ (1-K3)+K3 ZVy, (4.2) 
Vx Vx,e-™ (1-K3)+K3 ZVx, 


t+} 


® update-track-head: The position and velocity computed in center tangent and 
the new variance values are used as before to update the track-head. The state 
consists of the position latitude and longitude and the corresponding velocity 
components, while the covariance consists of the three steady-state values of the 
filtered variance. 

@® update-track-tail: This is modified in the same way as update-track-head to 
account for the new form of the state and covariance. 

The result is a close approximation to the variance values obtained from the 
full implementation and a position that represents the result of some combination of 
predicted position and contact position that was not specifically identified and operated 
on. The distribution of estimated filtered states closely approximates that of actually 


computed states. The specifics of this comparison are detailed below. 
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B. ACCURACY OF THE MODIFIED TRACKER 


The assumptions underlying this new procedure introduce deviations from the 
full implementations. As discussed above, the error in the positioning of the SPA 
about the target’s true position error is maximum for Ka approximately .5-.65. The 
mean value of this error, assuming the worst case of a non-maneuvering target, can 
be calculated using Equation (3.3) to find distance from the old filtered position to the 
new observation position. Since the observations are centered on the true target 
position, the mean of this value can be calculated using the true target position. The 
results are shown in Appendix E for favorable and unfavorable conditions. From 
these results it can be seen that the error is excessive for targets at high speed if the 
time between detections is very long. This quickly drops off however, and for more 
reasonable time intervals of less than one hour, the relative error, even for high speed 
targets, is negligible. 

Assumption three must also be valid for the approximation to be reasonable. 
The validity of the first part of the assumption 1s shown above. This requires that the 
variance of the distribution of filtered positions be equal to the filtered variance. 
Appendix F presents a comparison of filtered position distributions and filtered AOU 
sizes for various target speeds, computed using the model in the Appendix G. The 
data in Appendix F is the result of modelling 48 hours worth of tracking data or no 
more than 300 observations. The simulated target conducted a random tour consisting 


of ten course changes executed at time intervals with an exponential distribution and a 


51 


mean of 4 hours. The target’s speed was either 3 knots, 10 knots or 30 knots. The 
standard ASSET IOU parameters 8=.25 and o=106'” were used. The results of 30 
different iterations were averaged to determine the mean values. 

The results indicate that for long time intervals between contacts, and speeds 
near the modeled one, the results are fairly good. Excessive deviations are exhibited 
when the target’s speed is far in excess of the modeled speed. Even in the case 
where the actual target speed was ten knots, the damping of the velocity is apparent in 
the filtered mean speed. When this is combined with the offset error described above, 
errors of up to sixteen percent can result. While this is a large error, the average 
value this error takes through the range of common time intervals between 
observations (.1 hours to 1 hour) and possible contact speeds, is closer to five 
percent. 

One source of this error can be seen by comparing (3.10) and (3.12). While the 
single Kalman gain Ka is used to compute the value of the variance, the actual 
distribution is computed using both Ka and KbJ. For the range of AOU sizes 
associated with submarines and MPA assets (at worst 25nm) leaving the SPA 
unadjusted amounts to an error in standard deviation of the estimated distribution of 
less than Inm. The error is typically conservative, increasing the SPA size by the 
error, and is largest when the time interval between contacts is small. This produces 
errors which, while a significant percentage of the AOU size, are actually quite small. 


A notable exception is the decrease in the standard deviation in situations where the 
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target speed is high, counteracting some of the IOU process’ tendency to 
underestimate the targets speed in this case. 

The failure of the estimate to account for the target’s velocity causes the 
resultant track to be erratic. This effect worsens as the distance traveled between 
observations gets small with respect to the AOU size. The actual filter exhibits 
consistant errors in this case, while the estimate zig-zags wildly around the true target 
track. Since the linearity of the track is not critical to any part of ASSET except 
possibly the correlator, the effect is minimal. For any single point the estimate 
produces, a track could exist which includes it. The random draw combines elements 
from all these hypothetical tracks into a single collection which has the proper 
pointwise relationship to the real track. Thus, while specific combinations of 
observation AOU size and observation interval result in poor performance, others 
result in improved performance. 

The performance of the modified filter is reasonably consistent with the original 
filter operation. A 3-5 nm error in a 25 nm AOU does not significantly impact the 
cueing or search efforts of prosecuting platforms and tracker performance is actually 
improved for cases where the target is transiting, moving at high speeds without 
maneuvering. The errors present in the modified tracker appear to have no negative 
impact on the operation of the rest of the simulation and while introducing moderate 
errors in tracking maneuvering targets, balances that with improved performance 


against non-maneuvering ones. Appendix G contains examples of the original filter’s 
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tracking ability and that of the estimate for a variety of sensor AOU sizes and mean 


interarrival times. 
C. LIMITING THE RANGE OF FILTER OPERATION 


An interesting effect apparent from Appendix F, is the limiting nature on the 
filtered position distribution of the contact AOU on the one hand and the predicted 
AOU on the other. This situation stems from the same effect that moves the filtered 
position from the observation position toward the predicted position. This can be 
seen directly by looking at the values of the Kalman gains as a function of range for a 
given time interval between observations. Examples of these relationships are shown 
in Appendix H. As can be seen from these graphs, if the observation AOU radius in 
nautical miles is limited to being less than 150 times the time interval in hours, the 
Ka value remains above the "flat" portion of the curves, approaching zero. 
Observations which arrive at the tracker such that the interval since the last contact is 
too small compared to the AOU size will be filtered using primarily the prediction 
information. The effect of the contact is largely ignored. This is the filter’s way of 
saying that the information value of the contact is not high enough to be considered of 
value in computing the filtered position. 

Such contacts can be ignored as there is little difference between the predicted 
position and variance and the resultant state and variance. A lower limit of 1nm on 
the observation is hardwired into the make-state and compute-mean-cov objects. If 


the time between contacts is too short, little is gained unless the AOU is tiny. The 
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Inm limit implies that the minimum time interval between contacts is around .01 
hours to allow the contact data to be assimilated into the track. If the contact makes 
no contribution to the track update, there is no reason to perform the filter operation. 
The effect of this can be seen in the data in Appendix F, where the predicted and 
filtered values are very close for small time intervals and large observation AOUs. 
Thus in dense contact environments, with contacts arriving at the fusion center at 
intervals of less than 30 minutes, any contact based on an AOU larger than 150 times 
the time since the last contact in hours, can be ignored and not filtered. This value of 
150 corresponds to a Kalman Gain limit of approximately .1. The exact value is 
rather arbitrary so long as it not too far into the portion of the curves that change 
quickly with variation in observed AOU size. This could be made a user input to 
limit the resolution of the tracker and speed execution time during scenario 
development and testing, and then reduced to provide more accurate track information 
when the data is actually collected. 

In similar fashion tracks which have not been updated in a long time, do not 
contain enough information to affect contacts that arrive with AOUs less than 5-10 
times the time interval since the last observation. If the value of Ka is very close to 
the "flat" portion approaching one, the predicted value is ignored in favor of the 
contact data. This allows similar saving of filtering operations in sparse track 
conditions with accurate sensors. If the time interval between contacts is expected to 


be greater than one hour and the contact AOUs are small, using the raw contact to 
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update the track is a reasonable approximation if the observation AOU size is smaller 
than five times the time between contacts in hours. 

When these two "gating" techniques are utilized, substantial numbers of filtering 
calls can be eliminated. This can be particularly true in scenarios covering large 
ocean areas with many possible targets. During the initial stages of the campaign, as 
search assets slowly get cues for very quiet enemy tactical platforms, detections will 
not be made until the enemy is at close range. This will result in observations with 
relatively small AOUSs coming in at long time intervals. As cues are dispatched and 
friendly forces converge on the enemy, even a modest number of platforms soon 
begin producing contact streams that arrive at rapidly decreasing intervals. Putting 
limits on when the filter is invoked will prevent the simulation from wasting efforts 


computing a SPA that does not really use much of the contact information. 
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V. CONCLUSION 


From a discussion of the construction and use of the Kalman filter based 
Maneuvering Target Statistical Tracker it has been shown that the tracker 
implementation in ASSET is not currently matched well to the types of inputs it 
receives. This results in the calculation of many duplicate values, the performance of 
unnecessary matrix manipulations and the calculation of values in situations when 
little useful information is gained. Improvements to the tracker aimed at decreasing 
the computation time which were discussed include: 

® Use of the equations detailed in (3.3) through (3.12) to compute the actual 
filtered variance and position. 
@ Filter the variance of the contact distribution about the target prior to conducting 

a random draw, thus producing a position drawn from a distribution 

approximating that of the actual filtered positions variance. Since the 

observations play no part in calculating this variance, appropriate values for the 


sensor AOUs can be precomputed for a range of expected time intervals. 


@ Limit the situations where the filter is actually used to those which will result in 
a meaningful amount of information being extracted. 


® Using planar estimations of the latitude and longitude computations rather than 
perform the conversions on a "spherical earth". 


When taken together these modifications greatly reduce the complexity of the 


calculations required to update a track after an observation. They can also be put 
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under user control so the level of fidelity can be adjusted based on the requirements 
of the architecture and the stage in development of the architecture the analyst is at. 
Areas for future work in this area include: 


@ Formulating a faster way of estimating the steady state gains for use in a table 
of pre-calculated filtered variances. 


@ Analysis of the filtered position distributions of other types of trackers to 
determine if a better method of approximation can be determined, one which 
better accounts for target velocity. 


@ Formulating a simple, recursive estimate of the Kalman gains for on-the-fly 
calculation. An attempt at this by using limiting values as the time interval got 
small yielded very good steady state results, but divergence and chaotic 
behavior for small AOU sizes (<10nm radius) precluded its use iteratively. 


@ An analysis of the algorithms in the correlator to determine if all of its 
algorithms are necessary, or if a simple probability of correct/false/no 
correlation could be substituted. If appropriate, determine what those 
probabilities should be. 


@ Modelling how the AOU size should change with respect to range and bearing 
from the searcher. A strict linear bearing error vs range is a simple way of 
solving the problem of constant sensor AOU size. The real situation is not 
quite so simple. This requires computation of the beam sizes and beam 
distribution of a towed array and the effect of TMA techniques on determining 
the position course and speed of the target. This would allow for proper 
calculation and filtering of the sensor’s variance in velocity measurement as 
well. 

The last example opens up a plethora of possibilities for analysis of the parameters 
required for input to ASSET. A compendium of data on everything from 
communication network modelling, to aircraft maintenance data analysis, to platform 


on platform engagement simulation, and so on. Well researched values for the 


present and expected future interactions of platforms and their supporting 
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infrastructures need to be produced if ASSET and simulations like it are to be 
available to the average user. 

ASSET has the potential to be a very effective analysis tool. A great deal of 
work must be expended to provide the necessary data to allow its practical use. In its 
present form it requires so much expert knowledge on so diverse a collection of 
topics, it is doubtful one person could perform meaningful analysis without an 
inordinate amount of time spent in research. Continued evolution in the functionality 
of the simulation will provide increased applicability and accuracy of results. 
Increased streamlining of the computational algorithms will be necessary as this 
occurs to keep it running on a desktop workstation in a reasonable period of time. 
Striving for improvement in performance however, without providing a ready source 


of data will most probably preclude widespread use of the simulation. 
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APPENDIX A: CHARACTERIZATION AND ESTIMATION OF ¢ AND @. 


The following equations represent the elemets 1n the transition matnx, phi, and the IOU process noise, Q. 
Approximation equations which can be used to more quicly calculate the values are also given. 
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The graphs which follow show the behavior of these functions over various time intervals. The variable “time” 
is equal to the elapsed time from the last contact's time of observation to that of the present contact, 1n hours. 
The estimated values are computed using the first two significant terms in the series expansion of the exponential 
terms. These estimates are very accurate for interarrival times of two hours or less. 
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APPENDIX B: APPLICABLE SOURCE CODE 
This appendix contains portions of the ASSET source code which would be modified 
to implement the changes described above. The remarks indicate the appropriate 
portions which would be deleted or changed. This is not an all inclusive listing, but 


contains the portions most affected by the modifications. 


ACT-MTST Wed, Apr 25, 1980 1 


OBJECT ACT-MTST 


ACT-MTST 1s an Odject LISP implementation of MIST, the 
Manuevering Target Statistical Tracker. 


Ses S888 S22 SC SF SS SSF SSF SF SF SF SF SF SF S888 SSF 2 FFT SSSSSSSCFV SFT STF FFF S88 SSFVF FT SESE SE SESS ST SST TSF 2 


(defobfun (exist act-mtst) (init-list) 
(usual-exist init-list) 
(have ‘dtor (/ pi 180)) 
(nave ‘alpha 0.25) 
(nave ‘sigma (sqrt (/ 100 2))) : 10 knot target 


(defobfun (create-object-links act-mtst) (osis) 
(have ‘idbm (ask osis 1dbm)) 


(defobfun (start-new-track act-mtst) (uc track-id) 
(let* ((mean-cov-list (compute-mean-covariance uc) ) 
(mean (first mean-cov-list)) 
(var (second mean-cov-list)) 
(alt (obu-contact-altitude uc) ) 
(uc-spatial-dat& (obu-contact-spatial-data uc)) 
(lat (obue-data-field-lat uc-spatial-data) ) 
(Ing (obuedata-field-ing uc-spatial-data) ) 
(cte-list (list uc)) 
track state-cov 
) 
(setq state-cov (make-obu-state-field 
:ref-lat lat :ref-ing ing 
sstate mean :covariance var) ) 
(return-from start-new-track 
(make-obuctrack :id track-id 
smumber-of-contacts 1 
shead-state-covariance state-cov 
shead-contact-id (obu-contacte-id uc) 
shead-spatial-data uc-spatial-data 
:stail-state-covariance state-cov 
stail-contact-id (obu-contact-id uc) 
:tail-spatial-data uc-spatial-data 
saltitude alt 
:contacts ctc-list 
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(defobfun (compute-mean-covariance act-mtst) (uc) 
{let* ((uc-spatial-data (obu-contact-spatial-data uc) ) 
(type (obu-data-field-type uc-spatial-data) ) 
(lat (obu-data-field-lat uc-spatial-data) ) 
(lng {(obu-data-field-lng uc-spatial-data) ) 
(prg {(obu-data-field-pbrg uc-spatial-data) ) 
(maj {(obu-data-field-maj uc-spatial-data) ) 
(min (obu-data-field-min uc-spatial-data) ) 
(cse (obu-data-field-cse uc~-spatial-data) ) 
(cse-unc (obu-data-field-cse-unc uc-spatial-data) ) 
(spd (obu-data-field-spd uc-spatial-data) ) 
(spd-unc (obu-data-fleld-spd-unc uc-spatial-data) ) 
(mean (make-array 4 :initial-element 0)) 
(var (make-array ‘(4 4) :initial-element 0)) 
stationary smaj smin cosb sinb cosb2 sinb2 
) 
(setq stationary (/ (* sigma sigma) (* 2 alpha))) 
(setf (aref mean 2) (* spd (cos (* dtor (- 90 cse))))) 
(setf (aref mean 3) (* spd (sin (* dtor (= 90 cse))))) 
(setf (aref var 2 2) stationary) 
(setf£ (aref var 3 3) stationary) 
(setq smaj (max maj 1)) 
(setq smin (max min 1)) 
(setq cosb (cos (* dtor brg))) 
(setq sind (sin (* dtor brg))) 
(setq cosb2 (* cosb cosb)) 
(setq sinb2 (* sinb sinb)) 
(setq smaj2 (* 0.25 smaj smaj)) 
(setq smin2 (* 0.25 smin smin)) 
(setf (aref var 0 0) (+ (* smin2 cosb2) (* smaj2 sinb2))) 
(setf (aref var 0 1) (* (= smaj2 smin2) sinb cosb)) 
(setf (aref var 1 0) (aref var 0 1)) 
(setf (aref var 11) (+ (* smin2 sinb2) (* smaj2 cosb2))) 
(return-from compute-mean-covariance (list mean var)) 


(defobfun (update-track actemtst) (uc track) 

(let* ((uc-spatial-data (obu-contact-spatial-data uc) ) 
(id (obu-contact-id uc) ) 
(type (obu-data-field-type uc-spatial-data) ) 
(lat (obu-data-field-lat uc-spatial-data) ) 
(lng (obu-data-field-lng uc-spatial-data) ) 
(Drg (obu-data-field-brg uc-spatial-data) ) 
(maj (oObu-data-field-maj} uc-spatial-data) ) 
(min (obu-data-field-min uc-spatial-data) ) 
(cse (obu-data-field-cse uc-spatial-data) ) 
{cse-unc (obu-data-field-cse-unc uc-spatial-data) ) 
(spd (obu-data-field-spd uc-spatial-data) ) 
(spd-unc (obu-data-field-spd-unc uc-spatial -data) ) 
(Oobs-time (obu-data-field-obs-time uc-spatial-data) ) 
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(track-spatial-data (obu-track-head-spatial-data track) ) 
(tail-spatial-data (obu-track-tail-spatial-data track) ) 
(last-update-time (obu-data-field-obs-time track-spatial-data) ) 
(tail-update-time (obu-data-field-obs-time tail-spatial-data )) 
(last-state-cov (obu-track-head-state-covariance track) ) 
(tail-state-cov (obu-track-tail-state-covariance track) ) 
(last-ref-lat (obu-state-field-ref-lat last-state-cov) ) 
(last-ref-lng (obu-state-field-ref-lng last-state-cov) ) 
State-cov-old state-cov-ctc state-cov-new delta-time 
contacts time-old mean-var-list pos ctce ctc-sd count 

) 


(setq delta-time (- obs-time last-update-time) ) 


(if (>= delta-time 0) 
(progn 
(setq state-cov-old (change-tangent lat lng 
(Loumotion delta-time last-state-cov) )) 
(setq state-cov-ctc (make-state type lat lng brg maj min cse cse-unc 
spd spd-unc) ) 
(setq state-cov-new (center-tangent (filter state-cov-old 
state-cov-ctc) )) 
(setq track (update-track-head track uc state-cov-new 
uc-spatial-data id)) 
(setq track (update-track-tail track) ) 
) 
(progn 
(if (or (> obs-time tail-update-time) 
(< (obu-track-number-of-contacts track) 
(ask ldbm contacts-from-head-to-tail))) 
(progn 
(setq contacts (append (obu-track-contacts track) (list uc))) 
(setq contacts (sort contacts #‘' (lambda (x y) 
(< (obu-data-field-obs-time (obu-contact-spatial-data x) ) 
(obu-data-field-obs-time (obu-contact-spatial-data 
Ya) 
(if (>= obs-time tail-update-time) 
(progn 
(setq state-cov-old tail-state-cov ) 
(setq time-old tail-update-time) 
) 
(progn 
(setq mean-var-list (compute-mean-covariance (first 
contacts) )) 
(setq lat (obu-data-field-lat (obu-contact-spatial-data 
(first contacts)))) 
(setq lng (obu~data-field-lng (obu-contact-spatial-data 
(first contacts)))) 
(setq state-cov-old (make-obu-state-field 
:ref~lat lat 
:ref-lng lng 
:state (first mean-var-list) 
:covariance (second mean-var-list))) 
(setq time-old (obu-data-field-obs-time 


(obu-contact-spatial-data (first 
contacts) ))) 
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t)) 


cse-unc 


(last contacts)) ) 


conzacts))))) 


) 


) 


) 


) 


) 


Wed, Apr 25, 1990 


(setq count (min (ask ldbm contacts-from-head-to-tail) 
(- (list-length contacts) 1))/J 


(dotimes 
(setq 


(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setg 


(i count 


pos (+ (obu-track-number-of-contacts track) 


ctc 
ctc-sd 
type 
obs-time 
lat 

ling 

brg 

ma} 


(setq min 


(setq 
(setq 


cse 
cse-unc 


(setq spd 


(setq 
(setq 
(setq 


(setq 


spd-unc 


nil) 


(nth pos contacts) ) 
(obu-contact-spatial-data ctc)) 
(obu-data-field-type 


(obu-data-field-obs-time 


(obu-data-field-lat 
(obu-data-field-lng 
(obu-data-field-brg 
(obu-data-field-maj 
(obu-data-field-min 
(obu-data-field-cse 


(obu-data-field-cse-unc 


(obu-data-field-spd 


(obu-data-field-spd-unc 
Gelta-ctime (- obs-time time-old) ) 


ctc-sd) ) 
ctc-sq) ) 
ctc-sd) ) 
ctc-sqd) ) 
ctc~-sd) ) 
ctce-sd) ) 
ctc-sq) ) 
ctc-sd) ) 
ctc-sq) ) 
ctc-sd) ) 
ctc-sd) ) 


state-cov-old (change-tangent lat lng 
(ioumotion delta-time state-cov-old) )) 


state-cov-ctc (make-state type lat lng brg maj min cse 


spd spd-unc) ) 


(setq state-cov-old (center-tangent 
(filter state-cov-old state-cov-ctc))) 
(setq time-old obs-time) 


) 


(setq track (update-track-head track uc state-cov-old 
(obu-contact-spatial-data (first 


(= count) 


(obu-contacteid (first (last 


(setq track (update-track-tail track) ) 


(return-from update-track track) 


(defobfun (update-track-head act-mtst) 
(1et ((contacts (append (obu-track-contacts track) 


) 


(setq contacts (sort contacts #' (lambda (x y) 


y)))))) 


(setf (obu-track-number-of-contacts track) 
1)) 
(setf (obu-track-head-state-covariance track) 


track) 


(< (obu-data-field-obs-time (obu-contact-spatial-data x)) 


(track uc state-cov spatial-data id) 
(list uc))) 


(obu-data-field-obs-time (obu-contact-spatial-data 


(setf (obu-track-head-spatial-data 
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sState-cov 


) 


track) spatial-data) 


(+ (obu-track-number-of-contacts 


i 
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(setf£ (obu-track-head-contact-id track) id ) 
(setf (obu-track-contacts track) contacts ) 
(return-from update-track~-head track) 


ew SO OSS OS SS OSS SSS SOS S SESS S28 S8 SFT S2 SS SS VS VF SSS SSS SS SS 2 SS 2S SSF 2 SS SS eee ee ees |S 


(defobfun (update-track-tail act-mtst) (track) 
(let* ((nctcs (obu-track-number-of-contacts track) ) 
(nead-sc (obu-track-head-state-covariance track) ) 


(head-sd (obu-track~head~-spatial-data track) ) 
(head-id (obu-track~head-contact-1id track) ) 
(tail-sc (obu-track-tail-state-covariance track) ) 
(tail-sd (obu-track-tail-spatial-data track) ) 
(tail-id (obu-track-tail-contact-id track) ) 


(tail-update-time (obu-data-field-obs-time tail-sd)) 
(contacts (obu-track-contacts track) ) 
type obs-time lat lng brg maj min cse cse-unc spd spd-unc 
mean-var-list new-tail pos-new-tail new-tail-sd delta-time 
state-cov-old state-cov-ctc state-cov-new 
) 
(setq contacts (sort contacts #' (lambda (x y) 
(< (obu-data-field-obs-time (abu-contact-spatial-data x)) 
(obu-data-field-obs-time (obu-contact-spatial-data 


y))))) 
) 
(1£ (<= netcs 1) 
(progn (setf (obu-track-taill-state-covariance track) head-sc) 
(setf (obu-track-tail-spatial-data track) head-sd) 
(setf (obu-track-tail-contact-id track) head~-id) 
(return-from update-track-tail track) 
) 
) 
(if (<= netcs (ask ldbm contacts~-from-head-to-tail)) 
(if (= tail-id (obu-contact-id (first contacts) )) 
(return-from update-track-tail track) 
(progn 
(setq mean-var-list (compute-mean-covariance (first contacts) }) 
(setq lat (obu-data-field-lat (obu-contact-spatial-data (first 
contacts)))) 
(setq lng (obu-data-field-lng (abu-contact-spatial-data (first 
contacts)))) 
(setf (obu-track-tail-state-covariance track) (make-obu-state-field 
:ref-lat lat 
sref-lng lng 
sstate (first mean-var-list) 
scovariance (second mean-var-list))) 
(sett (obu-track-tail-spatial-data track) (obu-contact-spatial-data 
(first contacts) )) 
(setf (obu-track-tail-contact-id track) (obu-contact-id (first 
contacts) )) 


(return-from update-track-tail track) 
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(setq pos-new-tail (- (obu-track-number-of-contacts track) 

(ask ldbm contacts-from-head-to-tail))) 
(setq new-tail (nth pos-new-tail (obu-track-contacts track) )) 
(setq new-tail-sd (obu-contact-spatial-data new-tail)) 


(setq type (obu-data-field-type new-tail-sd) ) 
(setq obs-time (obu-data-field-obs-time new-tail-sd) ) 
(setq lat (obu-data-field-lat new-tail-sd) ) 
(setq ing (obu-data-field-lng new-tail-sd) ) 
(setq brg (obu-data-field-brg new-tail~sd) ) 
(setq ma} (obu-data-field-maj new-tail-sd) ) 
(setq min (obu-data-field-min new-tail-sd) ) 
(setq cse (obu-data-field-cse new-tail-sd) ) 
(setq cse-unc (obu-data-field-cse-unc new-tail-sd)) 
(setq spd (obu-data-field-spd new-tail-sd)) 


(setq spd-unc (obu-data-field-spd-unc new-tail-sd)) 
(setq delta-time (- obs-time tail-update-time) ) 


(setq state-cov-old (change-tangent lat lng (ioumotion delta-time 


ai. —-s¢))) 


(setq state-cov-ctc (make-state type lat lng brg maj min cse cse-unc spd 


spd-unc) ) 

(setq State-cov-new (center-tangent (filter state-cov-old 
state-cov-ctc) ) ) 

(setf (obu-track-tail-state-covariance track) state-cov-new) 


(setf (obu-track-tail-spatial-data track) new-tail-sd ) 
(set£ (obu-track-tail-contact-id track) (obu-contact-id new-tail)) 
(prune-contacts-before-tail track) ; Get rid of contacts before tail. 


(return-from update-track-tail track) 


(deforfun (prune-contacts-before-tail act-mtst) (track) 

(tet ((netes (Oobu-track-number-of-contacts track) ) 
(contacts (obu-track-contacts track) ) 
(max-ctcs (ask ldbm contacts-from-head-to-tail) ) 
tail-position 

) 
(setq tail-position (= netcs max-ctcs) ) 
(1f (<= tail-position 0) 
(return-from prune-contacts-before-tail nil) 

) 
(dotimes (i tail-position nil) 

(setq contacts (remove-element contacts 0)) 
) 
(setf (obu-track-number-of-contacts track) max-ctcs) 
(setf (obu-track-contacts track) contacts) 
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(defobfun (ioumotion act-mtst) (time state-cov) 
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(let ((lat (obu-state-field-ref-lat state-cov) ) 
(lng (obu-state-field-ref-1ng state-cov) )} 
(mean (obu-state-field-state state-cov) ) 
(var (obu-state-field-covariance state-cov)) 
motion-matrices phi f new-mean new-var 
vartmp phitrn 

) 
(if (= time 0) 
(return-from ioumotion state-cov) 
) 
(setq motion-matrices (calc-mot-mat time) ) 
(setq phi (first motion-matrices) ) 
(setq f (second motion-matrices) ) 
(setq new-mean (matrix-multiply phi mean) ) 
(setq vartmp (matrix-multiply phi var )) 
(setq phitrn (matrix-transpose phi)) 
(setq new-var (matrix-add (matrix-multiply vartmp phitrn) f)) 
(setq phi new-var) 
(setq phitrn (matrix-transpose new-var) ) 
(setq new-var (matrix-scalar-multiply 0.5 (matrix-add phi phitrn))) 
(return-from ioumotion (make-obu-state-field 
:ref-lat lat 
sref-lng ling 
:state new-mean 
:covariance new-var) 


S SSS SS SFP Ses SSS SS SSF SF SSF SS SF SF SF SSF SF SF SF 8S SEF SSO eTeF TFT SOF eFTeF SFT SeF ees ee wT SOFT ee es wees ee eee oe 


(defobfun (calc-mot-mat act-mtst) (time) 
(let* ((invalpha (/ 1 alpha)) 

(invalpha2 (* invalpha invalpha)) 

(expalphat (exponential (* -1 alpha time))) 

(exp2alphat (* expalphat expalphat) ) 

(f (make-array '(4 4) :initial-element 0)) 

(phi (make-array '(4 4) :initial-element 0)) 

) 
(if (> time 0) 

(progn 

(setf (aref f 0 0) (* invalpha2 (- time (* invalpha 
(- (* 2 (- 1 expalphat)) (* .5 (- 1 exp2alphat))))))) 
(sett (aref f 11) (aref £ 0 0)) 
(setf (aref £ 0 2) (* invalpha2 (+ 0.5 (- expalphat) 
(* .5 exp2alphat)))) 

(setf (aref £2 0) (aref £0 2)) 
(setf (aref f 1 3) (aref f 0 2)) 
(setf (aref £31) (aref £0 2)) 
(setf (aref f 2 2) (* 0.5 invalpha (- 1 exp2alphat))) 
(setf (aref £3 3) (aref f 2 2)) 


th Fh Fe eH OD 


) 
(progn 
(setf (aref 
(== (* 
(setf (aref 
(setf (aref 


0 0) (* invalpha2 (- (- time) (* invalpha 

(- expalphat 1)) (* .5 (= exp2alphat 1))))))) 
11) (aref £0 0)) 

0 2) (* dnvalpha2 (+ (- 0.5) expalphat 


> FR ND PH 


Th 
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(* -l1 .5 exp2alphat)))) 
(setf (aref 20) (aref f O 2)) 
(setf (aref 1 3) (aref £ 0 2)) 
(setf (aref 31) (aref £ 0 2)) 
(setf (aref f 2 2) (* 0.5 invalpha (- exp2alphat 1))) 
(setf (aref 3 3) (aref f 2 2)) 


rh Fe ke Fe rh 


) 
) 
(setq £ (matrix-scalar-multiply (* sigma sigma) f)) 
(setq phi (identity-matrix 4)) 
(setf (aref phi 0 2) (* invalpha (- 1 expalphat))) 
(setf (aref phi 1 3) (aref phi O 2)) 
(setf (aref phi 2 2) expalphat) 
(setf (aref phi 3 3) (aref phi 2 2)) 
(return-from calc-mot-mat (list phi f)) 


(defobfun (exponential act-mtst) (x) 
(if {(< x =350) 
(return-from exponential (exp -350)) 
) 
(12 (25x 5350) 
(return-from exponential (exp 350)) 
(return-from exponential (exp x)) 


(defobfun (make-state act-mtst) (type lat lng brg maj min cse cse-unc spd spd-unc) 
(let ((mean (make-array 4 sinitial-element 0)) 
(var (make-array '(4 4) :initial-element 0)) 
(stationary (/ (* sigma sigma) (* 2 alpha))) 
(smaj (max maj 1)) (smin (max min 1)) 
) 
(setf (aref mean 2) (* spd (cos (* dtor (- 90 cse))))) 
(setf (are= mean 3) (* spd (sin (* dtor (- 90 cse))))) 
(setq var (make-var stationary brg smaj smin)) 
(return-from make-state (make~-obu-state-field 
:ref-lat lat 
:refelng lng 
:State mean 
:covariance var)) 


5 —e ee SSS 2S 22 SS 2 2 22 2 2 2S OSS SESE SS2 OSS 2 @ 2 2 @ BS @ISS @ISSSP(SSPSSSSSSOSP(NAI*OSOSGSGOE OEE STS 22 == 


(defobfun (make-var act-mtst) (stationary brg smaj smin) 
(let* ((var (make-array '(4 4) :initial-element 0)) 
(cosb (cos (* dtor brg))) 
(sind (sin (* dtor brg))) 
(cosb2 (* cosb cosb)) 
(sinb2 (* sinb sinb)) 


2 


ACT-MTST 


(sma 42 
(smin2 
) 


(* 90.25 
(*0'.25 


(setf (aref var 
(setf (aref var 
(setf (aref var 
(setf (aref var 
(setf (aref var 
(setf (aref var 
(return-from make-var var) 


Wed, Apr 25, 1990 


smaj sma})) 

smin smin)) 

0 0) (+ (* smin2 cosb2) (* smaj2 sinb2))) 
O 1) (* (= sma 32 smin2) sinb cosb)) 

1 0) (aref var 0 1)) 

11) (+ (* smin2 sinb2) (* smaj2 cosb2))) 
2 2) stationary) 

3 3) stationary) 


° 


(defobfun (filter act-mtst) 


(state=-cov-sol state-cov-ctc) 


(let ((lat (obu-state-field-ref-lat state-cov-sol) ) 
(lng (Oobu-state-field=-ref-1ng state-cov-sol) ) 
(meanl (obu-state-field-state State-cov-sol) ) 
(varl (obu-state-field-covariance state-cov-sol) ) 
(mean2 (obu-state-field-state state-cov-ctc) ) 
(var2 (obu-state-field-covariance state-cov-ctc) ) 
(infmatl (make-array ' (4 4))) 

(infmat2 (make-array ‘(4 4))) 
(infvecl (make-array 4)) 
(infvec2 (make-array 4)) 
(tmpmat (make-array ‘(4 4))) 
(tmpvec (make-array 4)) 
mean var 


) 


(setq infmatl 
(setq infvecl 
(setq infmat2 
(setq infvec2 


(matrix-invert varl)) 
(matrix-multiply infmatl meanl) ) 
(matrix-invert var2)) 
(matrix-multiply infmat2 meanz2) ) 


(setq tmpmat (matrix-add infmatl infmat2) ) 
(setq var (matrix-invert tmpmat) ) 

(setq tmpvec (matrix-add infvecl infvec2) ) 
(setq mean (matrix-multiply var tmpvec)) 
(return-from filter (make-obu-state-field 


) 


:ref-lat lat 


ref-lng lng 
State mean 
covariance var) 


2 — eo Se Se 2 2 2 2 2S S22 SS SSS 2 S22 SSF SF SF SSF SF SST SST SSSSSSSSFS SSBF ST SFF SS SGeS282 Se See eee 
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(defobfun (change-tangent act-mtst) (newlat newlng state-cov) 


(let ((oldlat 
(oldlng 
(mean 
(var 
(rot 
(FOCErn 
(newtmp 


(obu-state-field-ref-lat State-cov) ) 
(Obu-state-field-ref-lng State-cov) ) 
(obu-state-fleld-state state-cov) ) 


(obu-state-field-covariance state-cov) ) 
(make-array ‘(4 4) :initial-element 0)) 
(make-array '(4 4))) 
(make-array ‘(4 4))) 
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(new-mean (make-array 4 )) 
(new-var (make-array '(4 4))) 
latinglist rngbrglist cenlat cenlng deltabrg cosb sinb 


rngl 
) 
(setq 


(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 


(setf 
(setf 


(sett 
(setf 


(setf 
(setf 
(setf 
(sett 
(setf 
(sett 
(setf 
(setf 


(setq 
(setq 
(setq 


(setq 
(setq 
(setq 


rng2 rng3 brgl brg2 brg3 


rngl (sqrt (+ (* (aref mean 0) (aref mean 0)) 

(* (aref mean 1) (aref mean 1))}})) 
brgl (arctan (aref mean 1) (aref mean 0))} 
latinglist (getlatlng oldlat olding rngl brgl 0)) 
cenlat (first latlnglist)) 
cening (second latinglist) ) 
brgl (second (getrngbrg oldlat olding cenlat cenlng 0))) 
brg2 (second (getrngbrg cenlat cenlng oldlat olding 0))) 
rngbrglist (getrngbrg newlat newlng cenlat cenlng 0)) 
rng3 (first rngbrglist) ) 
brg3 (second rngbrglist) } 
brg4 (second (getrngbrg cenlat cenlng newlat newlng 0))) 
deltabrg (mod (+ (= brg2 brgl) (- brg3 brg4) 360) 360)) 
cosb (cos (* dtor deltabrg) )} 
sind (sin (* dtor deltabrg})) 


Move the mean position: 
(aref new-mean 0) (* rng3 (cos (* dtor (= 90 brg3))))) 
(aref new-mean 1) (* rng3 (sin (* dtor (= 90 Brg3))))) 


Rotate the velocity: 


10 


(aref new-mean 2) (+ (* cosb (aref mean 2)) (* sinbd (aref mean 3)))) 
(aref new-mean 3) (- (* cosb (aref mean 3)) (* sinb (aref mean 2)))) 


Rotate the covariance: 
(aref rot 0) cosb) 
(aref rot 1) cosb) 
(aref rot 2) cosb) 
{aref rot 3) cosb) 
(aref rot 1) sinb) 
(aref rot 3) sinb) 
(aref rot 0) (= sinb)) 
(aref rot 2) (= sinbd)) 


WrNOWNF © 


newtmp (matrix-multiply rot var)) 
rottrn (matrix-transpose rot) ) 
new-var (matrix-multiply newtmp rottrn) } 


rot new-var) 
rotirn (matrix-transpose new-var) } 
new-var (matrix-scalar-multiply 0.5 (matrix-add rot rottrn))) 


(return-from change-tangent (make-obu-state-field 


tref-lat newlat 
:ref-lng newlng 
:state new-mean 


:covariance new-var) 


74 


ACT-MTST Wed, Apr 25, 1980 11 


(defobfun (center-tangent act-mtst) (state-cov) 
(let ((lat (obu-state-field-ref-lat state-cov) ) 

(lng (obu-state-field-ref-lng state-cov) ) 
(mean (obu-state-field-state state-cov) ) 

(var (obu-state-field-covariance state-cov) ) 

(rot (make-earray ‘(4 4) :initial-element 0)) 
(rottrn (make-array ‘(4 4))) 

(newump (make-array ‘(4 4))) 

(new-mean (make-array 4 )) 

(new-var (make-array ‘(4 4))) 

latlinglist rngbrglist newlat newlng newrng newbrg 
rngl brgl cosb sinb 
) 
(setq rngl (sqrt (+ (* (aref mean 0) (aref mean 0)) 

(* (aref mean 1) (aref mean 1))))) 
(if (< rngl 0.01) 
(return-from center-tangent state-cov) 

) 
(setq brgl (arctan (aref mean 1) (aref mean 0))) 
(setq latlnglist (getlatlng lat lng rngl brgl 0)) 
(setq newlat (first latlnglist)) 
(setq newlng (second latlnglist)) 
(setq rngbrglist (getrngbrg newlat newlng lat lng 0)) 
(setq newrng (first rngbrglist)) 
(setq newbrg (mod (+ (second rngbrglist) 180) 360)) 
(setq deltabrg (- newbrg brgl)) 
(setq cosb (cos (* dtor deltabrg))) 
(setq sinb (sin (* dtor deltabdrg))) 


Move the mean position: 
(setf (aref new-mean 0) 0) 
(setf (aref new-mean 1) 0) 


: Rotate the velocity: 
(setf (aref new-mean 2) (+ (* cosb (aref mean 2)) (* sinb (aref mean 3)))) 
(setf (aref new-mean 3) (- (* cosb (aref mean 3)) (* sinb (aref mean 2)))) 
Rotate the covariance: 
(setf (aref rot 0 0) coab) 
(setf (aref rot 11) coab) 
(setf (aref rot 2 2) cosb) 
(setf (aref rot 3 3) cosb) 
(setf (aref rot 0 1) sinb) 
(setf (aref rot 2 3) sinb) 
(setf (aref rot 1 0) (- sinb)) 
(setf (aref rot 3 2) {(- sinb)) 


(setq newtmp (matrix-multiply rot var)) 
(setq rottrn (matrix-transpose rot) ) 
(setq new-var (matrix-multiply newtmp rottrn)) 


(setq rot new-var) 
(setq rottrn (matrix-transpose new-var) ) 
(setq new-var (matrix-scalar-~-multiply C.5 (matrix-add rot rottrn))) 
(return-from center-tangent (make-obu-state-field 
:ref-lat newlat 
:ref-lng newlng 
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‘state hewmen 
:covariance new-var) 


(defopfun (arctan act-mtst) (y x) 


(let (a 
) 
(1f 


) 
(if 


(if 


) 
(if 


) 


S See eee eae es 


(and (= y 0) (# x 0)) 
(return-from arctan 0) 


(= y 0) 

(progn (if (< x 0) 
(return-from arctan 270) 
(return-fram arctan 90) 


) 


(< y 0) 
(setq a (+ (/ (atan (/ x y)) dtor) 180)) 
(setq a (/ (atan (/ x y)) ator) ) 
(> a 0) 


(return-fram arctan a) 
(return-fram arctan (+ a 360)) 


This procedure converts a covariance matrix (cov) into an 
ellipse. The ellipse is returned in the form: 


( Orientation , Semi-major axis , Semi-minor axis }. 


SOS Bae SOS Ge SS PS SESS OSS E88 SSF@PWSISS FSO OSS O88 SATS See eas 


(defobfun (cov-to-ellipse act-mtst) (cov) 


(let ((a 
(Dd 
(¢ 


(aref cov 0 0)) 
(aref cov 0 1)) 
(aref cov 1 1)) 


brg cosb sinb cosb2 sinb2 sincos2b smjr smnr 


) 
(if 


) 
(12 


) 
) 


(not (> (= (* ac) (* BD)) Q)) 
(return-from cov-to-ellipse (list 0 (expt 10 9) (expt 10 9))) 


(= Bb 0) 
(1f£ (not (> ac)) 
(return-from cov-to-ellipse (list 0 (* 2 (sqrt (abs c))) 
(* 2 (sqrt (abs a))))) 
(return-from cov-to-ellipse (list 90 (* 2 (sqrt (abs a))) 
(* 2 (sqrt (abs c))))) 


(1£ (> Bb 0) 
(4f (= a c) 


(return-fram cov-to-ellipse (list 45 (* 2 (sqrt (* a b) 


)) 
(* 2 (sqrt (- a bB))))) 


) 
{1f (< b 0) 
(if (= a c) 


) 

) 

(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 
(setq 


(return-fram cov-to-ellipse (list 135 (* 2 (sqrt (- a b))) 
(* 2 (sqrt (* a b))))) 


brg (/ (arctan (- ¢ a) (* 2 B)) 2)) 

sinb (sin (* dtor brg))) 

cosb (cos (* dtor brg))) 

sinb2 (* sinb sinb)) 

cosb2 (* cosb cosb)) 

sincos2b (* 2 Bb sinb cosb)) 

smjr (* 2 (sqrt (* (* @ sinb2) sincos2b (* ¢ cosb2))))) 

smnr (* 2 (sqrt (+ (* a cosb2) (* -1 sincos2b) (* ¢ sinb2))))) 


(return-from cov-to-ellipse (list brg smjr smnr)) 
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;lisp matrix manipulation routines 
74/28/89--s.c. lent 


; these routines operate on three(3) ‘standard type’ lisp objects which are 
subsets of lisp object types: 


1. scalars, i.e., lisp type ‘number' 
2. vectors, i.e., lisp type 'simple-vector’ with 


length > 1 
3. matrices, i.e., lisp type ‘array’ with 
rank = 2 and 


each dimension > 1 


Be te te Ye Ye Fe Se Ye We Ye 


all routines accept input, and produce output of these types 


unacceptable inputs produce ‘'NIL' as output. there are two(2) cases: 
1. "non-standard'’ inputs 
2. ‘standard' but mismatched inputs, e.g., 
the inverse of matrix(m x n), where m © n. 


te Se te ey Yo Ge %e 


;routines consist of: 
mx-mul t 
mx-add 
sc-mx-mult 
mx-trans 
mx-ident 
mx-inverse 


DUM hm WAH 


these are built using the following ‘sub-routines': 


matrixp 

. array-row-el 
array-col-el 
array~row-dim 
array-col-dim 
can-mx-mult 

. can-mx-add 
sqr-mx 
mx-factor 
mx-subst 


OwOdI nu AWAD! 


~ 


e Se Fe Fe Fe Fe Be Be Ce Bs eo We He Be Be Be Be Se Pe Fe Be 


e 


eeeeven ° ° ° ° eeeeene 
eoerveeeeveoeweveeeee 


(defun matrix-multiply (ab 
(if (not (can-mx-mult a b) 


(return-from matrix-multiply NIL) 


2 oo ete oeveeenvneoeseereeeevoeeeeve et eevee 
’ oorvrervoeeeosvseoweroevrevrwroererwrowoereoeH ee 


a ter Oe 


) 
(let ((¢) (ij-list)) 
(setq c (make-array (remove 1 (list (array-row-dim a) (array-col-dim b))) 
) 
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) 
(do ((indexl 0 (+ indexl 1))) 
((> indexl (- (array-col-dim b) 1))) 
(if (> (array-col-dim b) 1) (setq ij-list (list index1)) 
(setq ij-list ()) 
) 
(do ((index2 0 (+ index2 1))) 
((> index2 (- (array-row-dim a) 1))) 
(if (> (array-row-dim a) 1) (setq ij-list (cons index2 ij-list)) 
) 
(setf (apply #’aref c ij-list) 0) 
(do ((index3 0 (+ index3 1))) 
((> index3 (- (array-dimension b 0) 1))) 
(setf (apply #’aref c ij-list) (+ (apply #’aref ¢c ij-list) 
(* (array-col-el a (list index2 


index 3) ) 
(array-row-el b (list index3 


indexl)) 


) 
) 
(if (> (array-row-dim a) 1) (setq ij-list (cdr ij-list)) 
) 
) 
) 
(if (equal (array-dimensions c) () ) 
(aref c) 


(defun matrix-scalar-multiply (sce mx) 
(if (not (and (nmumberp sc) (matrixp mx))) 
(return-from matrix-scalar-multiply NIL) 
) 
(let ((ij-list) (mx-new (make-array (array-dimensions mx)))) 
(do ((indexl 0 (+ indexl 1))) 
((> indexl (- (array-col-dim mx) 1))) 
(if (> (array-col-dim mx) 1) (setq ij-list (list indexl)) 
(setq ij-list ()) 
) 
(do ((index2 0 (+ index2 1))) 
((> index2 (- (array-dimension mx 0) 1))) 
(if (> (array-dimension mx 0) 1) (setq ij-list (cons index2 ij-list)) 
) 
(setf (apply #’aref mx-new ij-list) (* (apply #’aref mx ij-list) sc)) 
(setq ij-list (cdr i4-list)) 


mx-new 


78 


Matrix functions Wed, May 1, 1988 


*@seeveeeveenvneevveeveeeeeveeeeevpeeeeveeeeeeeerseeveeeveeveeveeveeeveeeeeeeeeeseeoseneeveeveeeveeeeeeeneeee ee 
CEshstActeE sds isksisnciciclsicisibleis cislenelercicisieleloncielsiele cisivionsicicic cis ce eicleiale bee. s 


(defun matrix-transpose (mx) 
(1f (not (matrixp mx) ) 
(return-from matrix-transpose NIL) 
) 
(let ((ij-list) (mx-new (make-array (reverse (array-dimensions mx))))) 
(do ((indexl 0 (+ indexl 1))) 
((> indexl (- (array-col-dim mx) 1))) 
(1f (> (array-col-dim mx) 1) (setq ij-list (list indexl)) 
(setq ij-list ()) 
) 
(do ((index2 0 (+ index2 1))) 
((> index2 (- (array-dimension mx 0) 1))) 
(1£ (> (array-dimension mx 0) 1) (setq ij-list (cons index2 1j-list)) 
) 
(setf (apply @¢'aref mx-new (reverse ij-list)) (apply @'aref mx ij-list)) 
(setq ij-list (cdr ij-list)) 


eeseveeveeveeverer ee eevee eerevpeeeeeeeeeeeeveeeeeseeeseeseeeeeeseseeeseeeeseeeeee* se Goeovree ees 
(a 2 2 A 


(defun matrix-add (mxl mx2) 
(if (not (can-mx-add mxl mx2)) 
(return-from matrix-add NIL) 
) 
(let ((ij-list) (mx-new (make-array (array-dimensions mxl)))) 
(do ((indexl 0 (+ indexl 1))) 
((> indexl (- (array-col-dim mxl) 1))) 
(1f (> (array-col-dim mxl) 1) (setq ij-list (list indexl)) 
(setq ijelist ()) 
) 
(do ((index2 0 (+ index2 1))) 
((> index2 (- (array-dimension mxl 0) 1))) 
(if (> (array-dimension mxl 0) 1) (setq ij-list (cons index2 ij-list)) 
) 
(setf (apply @‘’aref mx-new ij-list) (+ (apply #‘aref mxl ij-list) 
(apply #*aref mx2 ij-list))) 
(setq ij-list (cdr ij-list)) 


eeeseese8s8? eeoseeeeseerss ee 
Par ar a a oe te 


se 


ee 
° 


eeoeeoee 
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(defun identity-matrix (m) 
(1f (not (and (integerp m) (> m 1))) 
(return-from identity-matrix NIL) 
) 
(let ((ij-list) (mx-new (make-array (list m m)))) 
(do ((indexl 0 (+ indexl 1))) 
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((> indexl (- m1))) 
(setq ij-list (list indexl)) , 
(do ((index2 0 (+ index2 1))) 
((> index2 (- m 1))) 
(setq ij-list (cons index2 i4j-list)) 
(setf (apply @'aref mx-new ij-list) (if (= indexl index2) 1 0)) 
(setq ij-list (cdr ij-list)) 


oeeveeveoceeeeeeeeeevee eee eeesevrerereeereeeeeere @eeeeseeeeseeoeeeeeeeee eee eeeeeseveeseeee eee @# 
eoeoeorvpeeveev ev eevewvrerevnevrereereeveereeeeerere ov eee ee e e e 


(defun matrix-invert (a) 
(if (not (sqr-mx a)) 
(return-from matrix-invert NIL) 
) 
(let* ((n (array-dimension a 0)) 
(temp (mx-factor a n)) 
(w (car temp) ) 
(t{pivot (cadr temp)) 


) 
(.£ (= 0 (caddr temp) ) 
(return-from matrix-invert NIL) 
) 
(m-subst wn ipivot) 


oweregeveveeveerervnererererevevneveeHeeoHweerHeeovn ea eoea eee eee eoeoevvvveev evr eveeveHerevrwvreeeevnevneHeeeee 


(defun matrixp (a) 
(if (arrayp a) 
(cond ((simple-vector-p a ) 
(4f (> (array-dimension a 0) 0) 
T 


NIL 
) 
) 
( (equal (array-rank a) 2) 
(4f (and (> (array-dimension a 0) 0) 
(> (array-dimension a1) 0) 


) 
é % 
NIL 
) 
) 
(T NIL) 
) 
NIL 
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(defun array-row-el (a i) 
(apply @’aref a (if (= (array-rank a) 1) 
(list (car i)) 


Pre ee ee ee ee ee ee ee 
rf ‘ ’ 


(defun array-col-el (a 1) 
(apply #’aref a (if (= (array-rank a) 1) 
(car 1) 


eseee*@eeoa*tPeaeewseweeaeeeeseepeeee*eeoeeee ee eeeetwoe eevee ese#eeseseeeeee*seseeeeeveeeeee*eteee eee 
i a ed 


(defun array-row-dim (a) 
(Lf (= (array-rank a) 1) 
Z 


(array-dimension a 0) 


(defun array-col-dim (a) 
(Lf (= (array-rank a) 1) 
a 


(array-dimension a (- (array-rank a) 1)) 


@e@oeoeee*@eeeostseoceeoeeveeeteewveeoeeeeeeeweeeeeeeeeoeeeee eeee*eeoeeveveeeeveeveet+e*aea ee toe se eevee eoeee ee 


(defun can-mx-mult (a b) 
(if (and (matrixp a) (matrixp b)) 
(1f (= (array-dimension a (- (array-rank a) 1)) 
(array-dimension b 0) 


) 
T 
NIL 
) 
NIL 


eoeeeet*@eoeee*e ee eee wae ene eaeeen eet eevese ee see ea ese see eseeeeeevneee*e eee eoneaee eee eee eet ee 6 
CORO EPO FOOEEOTOTFEOE EERO OECOEEEOC COED EDF OOo e TEL Oo EesoororsS 


(defun can-mx-add (a b) 
(if (and (matrixp a) (matrixp b)) 
(if (and (= (array-dimension a 0) 
(array-dimension b 0) 
) 
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(= (array-dimension a (= (array-rank a) 1)) 
(array-dimension b (- (array-rank b) 1)) 


) 
) 
T 
NIL 
) 
NIL 


(defun sqr-mx (a) 
({1f (matrixp a) 
(if (= (array-rank a) 2) 
(if (= (array-dimension a 0) (array-dimension a 1)) 


x 
NIL 
) 

) 

NIL 


ove , 2 ee ee ed 
oeorvree# st oer oer eoeowoneovreoewnewvrevnevreeeveevreeveaerowm eRe eoeeeoev eee wR evo eRe eE e 


(defun mx-factor (mx n) 
(let ((awikod) (colmax) (ratio) (rowmax) (temp) (istar) (iflag 1) 
(ipivot (make-array n)) 
(d (make-array n)) 
(w (make-array (array-dimensions mx))) 
) 
(if (<= n Q) 
(return-from mx-factor nil) 
) 


(do (( 
( 


QO (+ i 1)) 
owmax 0 Q) 


i 
E 
) 
((= i np) 
(setf (aref ipivot 1) i) 
(do “C(jeC ete. Je1))) 
((= j nd) 
(setf (aref wi j) (aref m i 4)) 
(setq rowmax (max rowmax (abs (aref wi 4)))) 
) 
(when (= rowmax 0) 
(setq iflag 0) 
(setq rowmax 1) 
) 
(setf (aref di) rowmax) 


eeeovrpevesoe eo 


(do ((k 0 (+ k 1))) 
Cie k (= nel) 
(setq colmax (abs (/ (aref w k k) (aref d k)))) 
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(setq istar k) 
(do (( 1 (+ m1) (+ £ 1))) 
((= £ n)) 
(setq awikod (abs (/ (aref wi k)(aref d 1)))) 
(when (> awikod colmax) 
(setq colmax awikod) 
(setq istar 1) 
) 


) 
(1£ (= colmax 0) 
(setq iflag 0) 
(progn 
(42 (> istar k) 
(progn 

(setq iflag (* -1 iflag)) 

(setq 1 (aref ipivor istar)) 

(setf (aref ipivot istar) (aref ipivot k)) 

(setf (aref ipivot k) 1) 

(setq temp (aref d istar)) 

(sett (aref d istar) (aref d k)) 

(secf (aref d k) temp) 

(do ((9 0 (+ 3 :1))) 
((= 9n)) 
(setq temp (aref w istar 4)) 
(setf (aref w istar 4) (aref wk 4)) 
(setf (aref wk 4) temp) 

) 

) 


) 
(do ((4 (* k 1p (* 2 1))) 
((2 i np) 
(setf (aref wik) (/ (aref wik) (aref wk k))) 
(setq. ratio (aref w i k)) 
(do (( 4 (> RB 1d 4 :1))) 
({* j np) 
(setf (aref wi 4) (~ (aref wi 4) (* ratio (aref wk 


J)))) 
) 
) 
) 
) 
) 
(14 (= (aref w (- nl)(=nil)) 0) 
(setq iflag 0) 
) 
(list w ipivot ifilag) 
) 
) 


eseeesevpee*ee fF @ 028 66 @ eoeeoee?e 
Serer rite cers @eeeoeeoeeeeeeeea eae eee eee ee ee 
Prerrrrerrerrrr rr rrr rere rere errr ere ee ere ee eee ee ee 


(defun mx-subst (wn ipivot) 


(let ( 
(D (make~array n :initial-element 0)) 
(ainv (Make~array (list nn))) 
(ip) 
) 


(do ((col 0 (* col 1))) 
((* col n)) 
(setf (aref b col) 1) 


eeoseeeven8eees * eecveeveeveveee ee eee ee 8 @ eoeeee 
° PMSA O Op OoOOORDOOOUO CO OOROCOOU OI re eas 


(setq ip (aref ipivot 0)) 
(setf (aref ainv 0 col) (aref b ip)) 
(do ((4 1 (* 2 1)) 
(sum 0 0) 
) 
((= 1 n)) 
(do ((9 0 (* 9 1))) 
((= 9 4)) 
(setq sum (+ sum (* (aref w i 3) (aref ainv 4 col)))) 
) 
(setq ip (aref ipivot i1)) 
(setf (aref ainv i col) (- (aref b ip) sum)) 
) 
(setf (aref ainv (- n1) col) (/ (aref ainv (~ n 1) col) (aref w (- n 


2) t=on Sy 
(do (( £ (e mn 2) (2 2 2)) 
( sum 0 0) 
) 
((# 1 -21)) 
(do (( 9 (# i 1) (+ 9 2D)) 
(f= jmp) 
(setq sum (+ sum (* (aref ainv j col) (aref wi 3)))) 


) 
(setf (aref ainv t col) (/ (- (aref ainv i col) sum) (aref wi :))) 


NOpDOdD TOC OODD OOS OODD ADHD OCR ar OOO RH Hee ea oo 
°. 7? * 
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Range-dearing functions Wed, Ape 25, 1980 1 


ae oe TRACK Utility Procedures: 


(defun getLatLng (latl ingl dist brg flag) 
iflag @ 0 great circle calc 
: = 1 rhumbline 
; Brg in degrees, 
(let* ((lat2 latl) (lng2 1lngl) (rng (rem (/ dist 60.1077) 360)) 
(cbrg (cosLisp brg)) (latit latl) (clatlt) (slatl) 
(srng nil) (crng nil) (ding nil) (temp nil) (clat2 nil) 


) 
(1f (* latl 90) (setq latit 89.9)) 
(af (@ latl -90) (setq latlt -89.9)) 
(setq clatl (cosLisp latlt)) 
(1f (= flag 0) > great circle calculation 
(progn (1f (> rng 180) 
(progn (setq rng (- 360 rng)) 
(setq brg (rem (* 180 Brg) 360)) 
(setq cbrg (cosLisp brg)) 
) 
) 
(12 (= clatl 0) 
(progn (1f (# (sinLisp brg) 1) 
(setq lat2 (- latl rng)) 
(setq lat2 (+ latl rng)) 


) 
(setq 1ng2 lngl) 
) 
(progn (setq slatl (sinlisp latlt)) 
(setq srng (sinLisp rng)) 
(setq crng (cosLisp rng)) 
(setq lat2 (- 90 
(acoslisp (fnAbs (+ (* clatl srng cbdrg) 
(* slatl crng)))))) 
(setq clat2 (cosLisp lat2)) 
(setq ding (acosLisp (fnAbs (/ (- crng (* slatl 
(sinLisp lat2))) 
(* clatl clat2))))) 
(42 (< (sinLisp brg) 0) 
(setq lng2 (fnCor (- 1lngl ding))) 
(setq lng2 (fnCor (+ 1lngl ding))) 
) 
) 
) 212 clati=d0 
) 3 progn rng > 180 
) :1f flag = gc 
(1f (* flag 1) srhumbline calc 
(progn (if (or (= Drg 90) 
(© brg 270)) 
(progn (setq lat2 latilt) 
(4£ (* brg 270) 
(setq lng2 (fnCor (- lngl 
(/ rng clatl)))) 
(setq lng2 (fnCor (+ lngl 


(/ rng clatl)))) 
) 
) send then construct brge=90 or brg=270 
(progn (setq lat2 (+ latit (* rng cbrg))) 
(12 (>= (abs lat2) 90) 
(progn (1f (> lat2 0) 
(setq lat2 90) 
(setq lat2 -90) 
) 
(setq lng2 0) 
) 
(progn (setq temp (log {(/ (tanLisp (- 45 (* 0.5 latlt))) 
(tanLisp (- 45 (* 0.5 lat2)))))) 
(setq lng2 (fnCor (+ lngl (* GLB. radToDeg temp 
(tanLisp brg))))) 
) 
) 
) ;end else construct brge90 or brge270 
) 
) 
) :; Lf Llagel 
(list lat2 1ng2) 
) let 
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(seftun getRngBrg (latl lngl lat2 lng2 flag) 
: flag = 0 great circle 
: = 1 rhumbline 
(let* ((rng 0) (org 0) (dlat (abs(- latl lat2))) 
(ding (abs(fnCor(- 1lngl 1lng2)))) (x) (y) (latlt) (clat2 (cosLisp lat2 )) 
) 


(i2 (>= (+ dat ding) 0.0001) 
(progn (if (* flag 0) 
(progn (if (= ding 0) 
(progn (if (< latl lat2) 
(setq brg 0) 
(setq brg 180) 
) 
(setq rng diat) 


(progn (if (= (abs latl) 90) 
(setq latlt (* 89.9 
(/ latl (abs latl)))) 
(setq latlt latl) 
) 
(setq clatl (cosLisp latlt)) 
(setq slatl (sinLisp latl)) 
(setq slat2 (sinLisp lat2)) 
(setq rng (acosLisp (fnAbs(+ (* slatl slat2) 
(* clatl clat2 
(cosLisp ding)))))) 
(setq brg (acosLisp (fnAbs(/ (- slat2 
(* slatl (cosLlisp 


rng))) 
rng)))))) 


(* clatl (sinLisp 


(42 (or (>= lng2 -90) 
(<= lngl 90)) <cother than east or west 


(progn (if (< lng2 lngl) 
(setq brg (- 360 brg)) 


dateline crossing 


: cessary . 
OnE ie are (progn (if (and (> lng2 90) 
; (< lngl -90)) 
(setq brg (- 360 brg))) 
») 


») 
)) 31 Llaged 
ae Tee (= (abs latl) 90) 
r i or (= s la 
ee (= (abs lat2) 90)) ;North or South pole 
(progn (setq rng dlat) 
(i{£ (> lat2 latl) 
(setq brg 0) 
(setq brg 180)) 


) 
x (1 / (tanti (- 68 (* 0.5 latl))) 
wT aot ae facies (= 65 (* 0.5 lat2)))))) 
(setq x (* x GLB. radTobdeg) ) 
(setq y (f£nCor (- lng2 lngl))) 
(setq brg (fnArctan x y)) 
(4f (or (= Dbrg 90) 
(= brg 270)) 
(setq rng (* clat2 ding)) 
(setq rng (abs (/ diat (cosLisp brg)))) 
) 
)) 
)) 
)) 
(setq rng (* 60.1077 rng)) 
(list rng brg) 
) slet 
) 
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Detection Reporter Wed, Mar 21, 1990 


eneaeoeoeeeeeneavnvae 
oeeoerev eer ero eroaevneoe” 


(defobfun (exist detectionReporter) (init-list) 
(usual-exist init-:ist) 


(have ‘report-type "POS+VEL ) ; POS-pos only, POS+VEL-pos + vel 
(have ‘positional-unc 10 ) ; 0 om 

(have ‘speed-unc 3 ) = kts 

(have *course-unc 39 ) ; degrees 


(have ‘expRepdelay 5) 
(have ‘FAR 0.9) 


(nave ‘starzTrackp 1) 


(sezq dataNamelist (append dataNameList 
'(:reports 
(report-type positional-unc speed-unc course-unct 
expRepDelay FAR startTrackp)))) 
(sets daztavValuelis= (append dataValuerist 
“(: reports 
(,repor=-type ,positional-unc ,speed-unc , Course-uns 
,expRepdelay ,.FAR ,startTrackp)))) 
(sezq dacaty;peLis=: (append datatypeiis= 
‘(:reports (data data data data data data data) ))) 
(setq datatexctlist (append dazaTextList 
*(:repores (°REPCRT TYPE" “POSITION ERROR” “SPEED ERRSR® 
"COURSE ERROR” “Exp Processing Dela;* 
“False A2arm Rate” "Track Inztiazticon”™)))) 
(sezq datatTemp.ateText (append dataTemplaceText 
*(: reports 
((°POS or POS“VEL® "(POS , POS*VEL!”) 
("Nastacas Mzles* "(0 = 10003%) 


(*Knots”® "(C = 33] >) 
(*Degrees® "(0 - 180)") 
(*"Hours” “(0 = LOO] 7) 
("FA per Day” "(0 - unlimited] *) 
("Single Contac: to Track” "0 - no, 1 - 
yes*))))) 
(lez* ((r5-designators ‘(reports-rb )) 
(ro-values (lisz dummyRSdata) ) 
(rd-names * (dutmynBdata) ) 


} 
(update-rb-lists rb-designators rb-values rb-names) 


ee J eoeeseeees ere eeeeesenee s seve eeeeteeesse ee 
wae siediciisilevedeieicne (el cies ‘steaelie eel sie. el eres 


(defoslun (setForStart detectionReporter) () 
(USSa-=SE@=5 S25 24-7) 
(nave ‘comsuteaFaAR FAR) 
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Detection Reporter Wed, Maer 21, 1980 


ee¢@e se 
SEB BEES E SOBRE ESES BESS AS OOOO CS HOCH DAC OOEOCHEODO OOH OOOO OO ON OOD pop oOo oon eee ete aoe 


eeseeeseeees 
opvpeoeoeoeeoed 


(defobfun (getReportDelay detectionReporter) (ésoptional targetType) 
(return-from getReportdelay 
(exponentialOraw expRepDelay) 


eevee eee 
opveoeeeawee 


(defobfun (make-obu-report detectionReporter) (time target-parms éoptional HFDFp) 


(let* ((cse (platformStateParms-cse target-parms) ) 
(sp2 (platfermScateParms-spd target-parms) ) 
(lat (pratformStateParms-lat target-parms)) 
(lag (platformStateParms-lng target-parms) ) 
(hs (platformStateParms-hgt target-parms) ) 
(targetName (platforrStaterarms-platformName target-pacms) ) 
(cat (make-obu-category :pinnedp nil 
slockedp nil 


:Cluster-count 0 
sassoc-count 93 o-0)) 
(sensor-field (make-obu-sensor 
:name (princ-to-string name) 
:Single-repost-to-trackp (not (equalp startTrackp 0)))) 
spartial-data-field 
) 


(if (equalp report-type 'POS+VEL) 
(setg cse (mod (+ cse (= (random (* 2.0 course-unc)) ccurse-unc) ) 
3€9) 
spd (max 0 (+ spd (- (random (* 2.0 speed-unc)) speed-unc)) 


(setyZ cse nil 
spd nil) 
) 
(multiple-value-setq (lat lng) (contactPositionDraw iat ling 
positional-unc) ) 


(setq spatial-data-field (make-obu-data-field 
type report -type 
:obs-time time 
slat lat 
zing ing 
:brg 0 
:ma4 positional-unc 
smin positional-unc 
:cse cse 
:¢se-unc course-unc 
: spd 
:spdeunc speed-unc 


) 
) 
(return-from make-obu-report (make-obu-contact 
ape | nil 
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Detection Reporter Wed, Mar 21, 1990 3 


sreceipt-time time 
:track-association targetName 

: Sensor sensor-field 
:categorization cat 

:Spatial-data spatial-data-field 
saltitude hgt 

: HFDFp HFDFp 


e@eeeeoee ee ee @eeseteeevee eee eve 
oo e?@ese oe ervovevreoev evo ee weve 


eve ° eet eoee oeeeee 
CCC oe hoo COO OE OOOOH EOHEFE OOH HOP OCHO H OO OOF ETO OCHERD Prereses 


° 
eoeoreeoevev ev ee 


(defobfun (make-false-obu-report detectionReporter) (lat lng uncertainty ctime 


targetType) 
(let ((falseTargetState (getFalsetargetState lat lng uncertainty targetType))) 


(make-obu-report ctime faiseTargetState) 
) 
) 


eeetsrre@ree eevee # ° eoeneee eee eevee eevee eevoveeveeeoseeveeevoevevee oe eevee ee 
@#eseopseeoeereoaservrveeaerevrwvrvrseseeeveowvneeaneeewvneve# oo e@# eo eve evovrevoorer eer vovoevrevrevreevresreaevrvr see oe wve 


(defobfun (getFalsetarsgezState detectionReporter) (lat lng uncertainty targetType) 
(muitiple-value-cind 
(lat lng) 
(contac=PositionDraw lat lng uncertainty) 
(return-from getFalseTarge:State 
(make-platformStateParms 
slat lat 
:lng lng 
“ngewe 
:cse (* (random 1.0) 369) 
:spd 10 
:platformName ‘FA 
:targetType tazge-Type 
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Structures 


oeoeoeoeoeeoeeteseeseseeeee* vee e*seevevnveeseoeveveeveeere Fees eveevseeveanccncveevseeeos 
oe 


Thu, bisy 10, 1980 


(defstruct coast-event 


ee er 
° 


time 
object 
procedure 
data 


updateList 


(defstruct coast-message 


eeenunvoeevae 
° 


send-time 
receipt-time 

type 

sender 

content 

size 
transmission-path 
transmission-count 


e®eseeseveveet+eseeeeeeeoeeseasane s 


[2 


(defstruct obu-contact 


) 


id 

receipt-time 
track-association 
sensor 
categorization 
spatial-data 
altitude 

EFDFp 


(defstruct obu-track 


) 


id 

number-of-contacts 
head-state-covariance 
head-contact-id 
head-spatial-data 
tail-state-covariance 
tail-contact-lid 
tail-spatial-data 
altitude 

(contacts nil :type list) 


(defstruct obu-data-field 


) 


type 
obs-time 
lat 

ing 

brg 

ma } 


min 
cse 
cse-unc 
spd 
spd-unc 


(defstruct obu-state-field 


) 


ref-lat 
ref-lng 
(state 


e@eetoeseoeeeseeeee+eeeve* eevee enna 


nil :type (array float 


; Time message sent. 
; Time message received at last node. 
> E.g. ‘Detect-msg 


*@eveeoeneeeveaeee ee 


4 )) 
)) 


(covariance nil :type (array float ‘(4 4) 


(defstruct obu-sensor 


) 


name 
single-report-to-trackp 


(defstruct obu-category 


*soesoe*@vove 


@eevevwosueee¢ @ e 


pinnedp 

lockedp 

assoc-count 

cluster-count 
ambiguity-list 
ambiguity-resolution 
ambiguity-resolution-time 


eeeees eefteeesee 
o7eeowevew vere e 


eee eeeetee 
[a 


eoaee 
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*#*@e5vreeeee*et#eoeeeeeeeeeeeee ee eva 
a a 


Probablitty functions Fri, Mar 23, 1980 


(defun contactPositionDraw (lat lng rng) 
; a draw from a circular normal -- rng is interpreted as a 2 sigma confidence 
interval 
(let ((randomRange (min (* (/ rng 2) (sqrt (* -2 (log (random 1.0))))) 
(* @ rng))) 
(randomAngle (randomCse) ) 
) 
(return-from contactPositionDraw 
(values-list 
(getLatLng lat lng randomRange randomAngle 
*greatCircle*))) 


(defun expandAOUpos (vel tau timeLate) 
(let ((ratio (/ timeLate tau)) 
(factor (* 2 pi vel vel tau tau))) 

(* factor (+ ratio -l (exp (* -1 ratio)))) 


eoeoeoeor ee? GO eee eeeerevreeveeoeveeoeneneoeoneoeoeeoneeeeene 


(defun expandAQUposVel (vel tau timeLate) 
(let* ((zatio (/ timeLate tau)) 
(factorl (* 2 pi vel vel tau tau)) 
(factor2 (* 2 (exp (* <1 ratio)))) 
(factor3 (* 0.5 (exp (* -2 ratio)))) 
) 
(* factorl (+ ratio factor2 -1.5 (* -1 factor3))) 


e®eeeeeowveeevreeevreeeeeeesr eevee eeeeeeeee ee eeesveeeeeeeeeeeeeveveeeeeeee*rereevreeeeervrseeeeane @ 
[2 


(defun getRandomPointInRegion (region) 
(let* ((regionParms (getRegionParms region) ) 
(minLat (region-parms-minLat regionParms) ) 
(maxLat (region-parms-maxLat regionParms) ) 
(mining (region-parms-minLng regionParms)) 
(maxLng (region-parms-maxLng regionParms)) 
(ckStuff (region-parms-checkpts regionParms) ) 
(thePoint nil) 
) 
(if (> mining maxing) (setq minlng (- minlng 360))) 
(do ((insidep nil)) 
(insidep) 
(setq thePoint (list 
(+ minlat (* (- maxlat minlat) (random 1.0))) 


(+ minlng (* (- maxing mining) (random 1.0))))) 
(setq insidep (checkIfInside thePoint ckStuff) ) 


) 
(return-from getRandomPointInRegion thePoint) 
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APPENDIX C: ITERATIONS TO ACHIEVE STEADY STATE GAIN VALUES 


The model in this appendix calculates the percent differece between succesive iterated values of Kalman gain K 
These values are tabulated for various ranges of interarrival time and sensor AOU size. 
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tN 
Nw 
th 


r, 
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APPENDIX D: KALMAN GAIN VERSUS INTERARRIVAL TIME 


This appendix shows the relationship between the Kalman gain and the interarrival time. This 
information can be used to select a range of time values to tabulate preprocessed variance data. 





















2 
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Vig = tl v2g = 0 V3BQ = 3 $2 = +-[1-emp[-a-A]] 3 = exp[-a-A]} 2 
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ne | a a) | [vt 42:$2-v2, 4 $2]°+v3, 4a! tri | ([ $3]°+v3, +q3 423 -[[ v2q +H2-v3, |- $3 492 | 
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2g Pe kta, = > Bp = B= > APPENDPRN(Kalman) =[A u klzs k2ar5 K2bo¢ k39¢ |e 
KDATA = READPRN(Kalman) = iene = KDATAD> KI = KDATAS> n = 0..laq KDATA 2 >| 
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l 
Kl 10.5 
8.001 0.0! 0.1 ! 10 


ume, 
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KDATA : : 
READPRN(KalmanS) | O> = KDATA > a = 0..last KDATA © > | 
Sensor AOU =S5nm 
K1 40.5 
8 O01 0.01 0.1 10 
time 
n 
KDATA = READPRN(Kalmani0 
(Kalmani0) ine = KDATA> KI = KDATA> n = 0..las§ KDATA™ > 
Sensor AOU =10nm 
K1 0.5 
SaaS eS 
8 ool 0.01 O01 10 
tume 
n 
KDATA = READPRN(Kalman25 
(Kalman25) ine = KDATA> Ki = KDATAZ> n = 0..las KDATA> | 
Sensor AOU =25nm 
K1,0.5 
01 01 , 10 
tume 


94 


KDATA - READPRN(KalmanS0 ) "erga i me 0 >! 
_ S ee J 


Sensor AOU =50nm 


01 10 100 
ume, 
KDATA = READPRN(Kalm100 
(Kalml00) me = KDATA2> ki = KDATAZ> n = 0..las{ KDATA >| 
Sensor AOU =100nm 
! 
K1,0.5 
Pon 10 100 
time 


APPENDIX E: OFFSET ERROR OF FILTER POSITION DISTRIBUTION 
This appendix provides indications of how far from the target’s true position the 
filtered position is displaced by the IOU process adjustment of velocity to account for 
the possibility of a manuever. This assumes the target is actually moving with 
constant course and speed. The first graph cover’s the entire range from inm to 


1000nm while succeding graphs show only the neighborhood of the maximum value 


Offset Error for 10 kt Target. 


Interarrival Time: 5 hours. 


12 


10 


Percent Error Compared to Fittered AOU. 
oO 


1 10 100 1000 
Sensor AOU 2-sigma radius (in nm). 
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Percent Error Compared to Fittered AOU. 


Percent Error Compared to Filtered AOU. 


Offset Error for 30 kt Target. 


Interarrival Time: 5 hours. 





Sensor AOU 2-sigma radius (in nm). 


Offset Error for 3 kt Target. 
Interarrival Time: 1 hour. 
0.314 


0.313 
0.312 
0.311 

0.31 
0.309 
0.308 
0.307 
0.306 
0.305 


10 105 11 #115 12 125 13 135 14 #145 15 
Sensor AOU 2-sigma radius (in nm). 
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Percent Error Compared to Filtered AOU. 


Percent Error Compared to Fittered AOU. 


Offset Error for 10 kt Target. 
Interarrival Time: 1 hour. 
1.045 


1.04 
1.035 
1.03 
1.025 
1.02 


1.015 
10 105 11 #115 12 125 13 135 14 145 15 


Sensor AOU 2-sigma radius (in nm). 


Offset Error for 30 kt Target. 
Interarrival Time: 1 hour. 
3.14 


3.13 
3.12 
3.11 
3.1 
3.09 
3.08 
3.07 
3.06 
3:03 


0 105 11 11.5 12 125 13 135 14 145 15 
Sensor AOU 2-sigma radius (in nm). 


Percent Error Compared to Filtered AOU. 


Percent Error Compared to Filtered AOU. 


Offset Error for 10 kt Target. 
Interarrival Time: 0.1 hour. 
0.0138 


0.01375 
0.0137 
0.01365 
0.0136 
0.01355 


0.0135 
1 105 1.1 1.15 1.2 125 13 #135 1.4 1.45 1.5 


Sensor AOU 2-sigma radius (in nm). 


Offset Error for 30 kt Target. 
Interarrival Time: 0.1 hour. 
0.0414 


0.0413 
0.0412 
0.0411 

0.041 
0.0409 
0.0408 
0.0407 
0.0406 
0.0405 


1 105 1.1 1.15 12 125 13 135 14 1.45 1.5 
Sensor AOU 2-sigma radius (in nm). 


APPENDIX F: ERROR BETWEEN FILTERED AND ESTIMATED SPA SIZE 
This appendix shows the difference in standard deviation between the actual 
filtered position distribution and the estimated distribution. The data tabulated shows 
more complete data for the cases modeled. This includes the observed AOU 
(Obs.AOU) and the predicted AOU (PRED AO); the mean value of the standard 
deviation calculated for 30 runs of 48 hours of 300 observations for the position 
predicted by the IOU model, the contacts position, and the filtered position; as well as 


the actual mean speed of the contact over the 30 runs and the mean filtered value. 


Time=4 Hours, Target S =10kts. 
Parameters: Alpha=.25, V= 10 kts. 





Time=1 Hour, Target Speed =3kts. 
Parameters: Aipha=.25, V= 10 kts. 





Time= : icy WLS Seta Okts. 





Time=1 Hour, Hour, Target pote =30kts. 





0 20 4) 60 80 100 120 140 160 


Sensor AOU Standard Deviation in nm 


emes) Poze we Esurnated Position 


Time=. | Hour, EU ete oe een 


Agpea) Pusan Gocasted Postion 
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Data for on 1OU process motion model beta=.25 and sigma based on 10 knot average speed 
Time interval between observations of 4 hours: 


Standard Deviation of error in nautical miles: 


Prediction Contact 


mean 
14.941 
20.041 
22.194 
33.145 
52.749 
92.721 
149.547 


21.938 
27.577 
27.27 
40.463 
57.402 
95.207 
146.205 


56.41 
59.391 
63.26 
78.533 
101.683 
133.696 
169.943 


mean 
11.957 
19.071 
23.148 
43.599 
81.301 
160.022 
245.976 


11.925 
19.357 
22.54 
43.917 
84.078 
154.341 
243.202 


13.4 
21.688 
26.691 
47.705 
89.644 

169.918 
233.559 


Fitter 

mean 
11.032 
16.716 
19.693 
32.113 
53.04 
93.447 
151.585 


11.197 
17.561 
19.703 
34.288 
94.921 
95.083 
147.258 


13.682 
22.186 
26.537 
47.533 
78.28 
121.254 
160.386 


SPA 

11.499 
17.18 
20.594 
35.113 
59.109 
103.544 
148.352 


11.499 
17.18 
20.594 
35.113 
59.109 
103.544 
148.352 


11.499 
17.18 
20.594 
35.113 
59.109 
103.544 
148.352 


%Error 
4.06122 
-2.70081 
-4.37506 
-8.54384 
-10.2675 
-9.75141 

2.179276 


-2.62632 
2.217695 

4.3265 
-2.34956 
-7.08522 
8.17141 
-0.73744 


18.98426 
29. 13853 
28.85792 
35.37151 

32.4333 
17.10384 
8.111788 


Time interval between observations of 2 hours: 


Standard Deviation of error in nautical miles: 


Prediction Contact 


mean 
7.08 
10.401 
16.528 
25.218 
43.418 
58.566 


9.704 
12.211 
19.499 
28.038 
43.761 
64.397 


mean 


4.859 
9.202 
20.859 
40.585 
78.444 
115.538 


4.805 
8.83 
21.68 
39.535 
78.011 
115.971 


Fitter 

mean 
4.547 
8.001 
15.156 
24.689 
43.725 
58.707 


4.583 
7.806 
16.535 
25.836 
42.761 
63.61 


SPA 
4.658 
8.471 

17.217 
28.317 
47.058 
64.879 


4.658 
8.471 
17.217 
28.317 
47.058 
64.879 


% Error 
-2.383 
-5.54834 
-11.9707 
-12.8121 
-7.08275 
-9.51309 


-1.61013 
-7.85031 

-3.9612 
8.76152 
-9.13128 
-1.95595 
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Obs. AOU Pred. AO mean 


25 
40 
50 
100 
200 
400 
600 


25 
40 
50 
100 
200 
400 
600 


25 
40 
50 
100 
200 
400 
600 


30.684 
34.216 
36.52 
47.234 
67.285 
108.788 
152.841 


30.684 
34.216 
36.52 
47.234 
67.285 
108.788 
152.841 


30.684 
34.216 
36.52 
47.234 
67.285 
108.788 
152.841 


Obs. AOU Pred. AO 


10 
20 
50 
100 
200 
300 


10 
20 
50 
100 
200 
300 


14.219 
17.145 
24.218 
33.625 
50.636 
67.655 


14.219 
17.145 
24.218 
33.625 
50.636 
67.655 


Speed in knots. 
Contact Fittered 
mean 
2.959 3.044 
2.975 3.145 
3.078 3.078 
3.02 2.635 
2.92 2.118 
2.992 1.881 
3.061 1.895 
10.133 7.482 
10.026 7.087 
9.997 6.81 
9.952 6.369 
9.943 5.999 
10.307 5.925 
10.076 5.958 
29.933 21.703 
29.931 20.036 
29.962 20.124 
30.019 18.783 
30.001 18.033 
29.908 17.635 
30.072 17.556 
Speed in knots. 
Contact Filtered 
mean mean 
3.041 3.374 
3.017 3.416 
3.023 2.904 
3.068 2.421 
3.065 2.144 
2.897 1.979 
10.066 8.383 
9.972 7.805 
9.964 7.004 
10.025 6.615 
9.958 6.384 
9.975 6.213 


22.67 5.196 5.242 4.658 1253757 
25.584 10.066 9.858 8.471 16.37351 
32.774 22.442 20.253 17.217 17.63373 
47.178 43.851 35.221 28.317 24.38111 
64.995 83.895 56.581 47.058 20.23673 
77.334 119.128 71.016 64.879 9.459147 


Time interval between observations of 1 hours: 
Standard Deviation of error in nautical miles: 


Prediction Contact __ Filter 
mean mean mean SPA %Error 

3.607 2.323 2.146 2.296 -6.5331 
5.215 4.445 3.803 4.159 -8.55975 
8.437 10.543 7.439 8.435 -11.8079 
11.984 19.878 11.43 13.716 -16.6667 
18.691 37.993 18.477 22.018 -16.0823 
25.46 58.34 25.401 29.352 -13.4608 
4.568 2.393 2.243 2.296 -2.30836 
6.018 4.461 3.846 4.159 -7.52585 
9.107 10.411 7.472 8.435 -11.4167 
13.545 19.994 12.392 13.716 -9.65296 
19.93 37.793 19.213 22.018 -12.7396 
25.759 54.988 25.225 29.352 -14.0604 
9.34 2.448 2.452 2.296 6.794425 
10.781 4.784 4.551 4.159 9.425343 
14.764 10.961 9.315 8.435 10.43272 
20.459 21.364 15.801 13.716 15.20122 
28.395 39.358 24.728 22.018 1230811 
36.583 57.941 33.572 29.352 14.37721 


Time interval between observations of .5 hours: 
Standard Deviation of error in nautical miles: 


Prediction Contact _ Filter 

mean mean mean SPA %Error 
1.03 0.488 0.474 0.485 -2.26804 
2.563 2.203 1.847 2.033 -9.14904 
3.622 4.095 3.014 3.05 -1.18033 
5.894 9.608 5.522 6.688 -17.4342 
8.469 18.734 8.256 10.542 -21.6847 
11.295 27.866 11.156 13.745 -18.8359 
1.533 0.487 0.475 0.485 -2.06186 
2.839 2.205 1.874 2.033 -7.82095 
4 4.182 3.131 3.5 -10.5429 
6.342 9.788 5.743 6.688 -14.1298 
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10 
20 
50 


200 
300 


14.219 
17.145 
24.218 
33.625 
50.636 
67.655 


Obs. AOU Pred. AO 


5 
10 
25 
50 

100 
150 


5 
10 
25 
50 

100 
150 


10 


25 


100 
150 


6.429 
8.15 
12 
16.684 
24.227 
31.113 


6.429 
8.15 
12 
16.684 
24.227 
31.113 


6.429 
8.15 
12 
16.684 
24.227 
31.113 


Obs. AOU Pred. AO 


1 
5 
10 
25 
50 
75 


1 
5 
10 
25 


2.162 
3.755 
5.499 
8.187 
11.759 
14.768 


2.162 
3.755 
5.199 
8.187 


30.061 
29.933 
30.047 
29.957 
29.938 
29.979 


24.285 
21.955 
20.111 
19.565 
18.814 
18.895 


Speed in knots. 
Contact Filtered 


mean mean 
3.079 3.489 
2.989 3.515 
2.96 3.119 
3.032 2.737 
2.984 2.329 
2.991 2.215 
9.988 8.922 
9.904 8.387 
10.002 La 
10.005 7.255 
10.091 7.009 
9.946 6.891 
30.02 26.103 
30.055 24.381 
29.998 22.467 
29.948 21.356 
29.953 20.68 
30.002 20.393 
Speed in knots. 
Contact Filtered 
mean mean 
3.046 3.178 
2.988 3.437 
3 3.298 
3.033 2.875 
3.028 2.58 
3.01 2.435 
10.011 9.607 
9.984 8.905 
9.944 8.461 
10.024 7.945 


9.101 18.558 8.7 10.542 -17.473 
11.972 27.983 11.653 13.745 -15.2201 
3.513 0.51 0.513 0.485 5.773196 
4.761 2.365 2.207 2.033 8.55878 
6.158 4.456 3.791 3.5 8.314286 
9.329 10.247 7.273 6.688 8.74701 
13.198 19.553 11.421 10.542 8.338076 
16.49 28.111 14.945 13.745 8.730447 
Time interval between observations of .1 hours: 
Standard Deviation of error in nautical miles: 
Prediction Contact _ Filter 
mean mean mean SPA %Error 
0.45 0.43 0.34 0.377 -9.81432 
1.155 1.921 1.059 1.263 -16.152 
1.689 3.728 1.621 1.995 -18.7469 
2.084 5.502 2.034 2.573 -20.9483 
0.495 0.429 0.348 0.377 -7.69231 
1.241 1.934 1.112 1.263 -11.9557 
1.761 3.748 1.665 1.995 -16.5414 
2.221 5.549 2.144 2.573 -16.6731 
0.778 0.466 0.428 0.377 13.52785 
1.721 2.005 1.368 1.263 8.313539 
2.408 3.832 2.104 1.995 5.463659 
2.947 5.618 2.685 2.573 4.352895 


Time interval between observations of .05 hours: 
Standard Deviation of error in nautical miles: 


Prediction Contact Filter 
mean mean mean SPA %Error 
0.322 0.4 0.274 0.312 -12.1795 
83 1.843 0.795 0.978 -18.7117 
0.359 0.41 0.291 0.312 -6.73077 
0.871 1.855 0.82 0.978 -16.1554 
0.534 0.442 0.37 0.312 18.58974 
1.176 1.908 1.045 0.978 6.850716 
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50 
15 


10 
25 


75 


11.759 
14.768 


2.162 
3.755 
5.199 
8.187 
11.759 
14.768 


Obs. AOU Pred. AO 


1 
5 
10 
15 


10 
15 


10 
15 


0.592 
1.504 
2.22 
2.78 


Obs. AOU Pred. AO 


1 
5 


1 
5 


Gg — 


0.407 
1.077 


0.407 
1.077 


0.407 
1.077 


10.05 7.625 
9.981 7.512 
30.01 28.592 
30.038 26.093 
30.015 24.905 
29.992 23.28 
29.652 22.652 
30.026 22.335 
Speed in knots. 
Contact Filtered 
mean mean 
3.002 3.192 
3.021 3.081 
3.003 2.847 
2.995 2.734 
10.003 9.568 
10.005 8.966 
10.026 8.711 
10.015 8.563 
29.986 28.414 
29.998 26.532 
30.037 25.88 
30.008 25.47 
Speed in knots. 
Contact Filtered 
mean mean 
2.971 3.138 
2.975 2.919 
9.996 9.549 
10.028 9.011 
29.997 28.449 
30.011 27.085 


APPENDIX G: MODELS OF ACTUAL ESTIMATED FILTER OPERATION 


This file simulates the basic operation of the ASSET tracker. A contact conducts a random tour at constant speed 
as a contact would in the simulation. Time interval between glimpses is an exponentially distributed variable. 
Statistics are collected to see how well the tracker tracks the target over the course of its path. The number of 
points plotted was adjusted to allow details to be seen. 

j = Otees runtime = 48 ghmpse = .1 A = -In(md(1)): glimpse m= 4° v = 10 tmp = 0 u = 10 mg - u 


tm, 4 = un, +4 al until runtime - tm; ,) N2 = lat(N)-1t = 0..N2) a = 0..N2 4 

Al = -In(md(1))*m A2 = -In(end(1))-m HAL =A. = - In(rnd(1))*m +42 A4 = -In(end(1))>mi +43 AS = -In(end(1))- mi 1A4 

Q6 = -In(md(1))*m +4547 = -In(end(1))-m +4648 = -In(end(1))-m +47A9 = -In(rnd(1))> mi +A8cl = 2+ Hemmd(1 )c2 = 2+ 4° mnd(1) 

c3 = 2¢H+md(1 )cd = 2+H+mnd(1 cS = 2-Wemd(1 )c6 = 2-Hemd(1 )c7? = 2°H°mnd(1 )cB = 2°Hemd(1) c9 = 2°H+md(L) clO = 2° H°mnd(1) 

te, = if um, <AL cl ,if tm, <A2,c2 if em, <Q3 3 vif im, <A4 ,c4 if tm, <O5 ,c5,if tm, <A6 ,6 ,if tm, <A7 7, if tm, <A8 ,<8 ,if um, <A9, 
rand|, = mnd(1) = rand2, = md(1) rand3, = md(1) rand4, = md(1) 


ong). > - =e, a : a ; =a . ; : : ‘ 
mg, - = °.t -2°In rand: cos 2° K-rand2, ymg, = Sad? In rand}, sin: 2 ys rand4,, | mse, = ee es 
rspd_, = v +{end(6)- 3) rorg, = 2° x: end(1) Yspd_, = vesini te, | Xspd, = Vecosi | tes Ydiay = 0 Xdisp =: 0 


L 


— 3s 
100 o 2 Bi r 1 
: - : = 25 =o 3 = — 1 = .25° ly, = tl v2, = 0 v3, = 03 > —* l-exp -a@ Li 
Yds, 4,’ | Yspd.-A, +¥dis, 7 3 oon | |) ac 0 oy i 
Xdis. 4,’ Xspd.-AL +Xdis,_ 2 - ij 
n a | SP n an n e 4 Oo f l rr e - 13 - r V4 | 
= Oe 1 = =" oe cca mc pt YS “O° eo oe pie Or * 2 ils 
Yspd vesin te 3 | OO | shir ioe bait. ar Qa! | L Pl et 
Xspd_ : vi cos ten : o | - 
lat, yg, +Ydis, q2, = aan he exp! -ar Ay, + -Svexp ar 2° A! 3, : Be ie exp! - ar 2° Ay 1 
Ing, xg, +Xdis, = 
yrspd, oa 
xrspd _fspd,* cos! f eae y | 
62. -v2. £462.17 Nl eater ae Manas i 
7 Vig 42° b2y°v2, + $2, | +v3, tal!) $3, *v3, ta3, +3 -.. v2, +O2,°v3, O35, H02, | 
; - ~2 Tre 2 one : 2 
VL, 42° O2y°v2, + G2, + v3, tql, tl a O35, (°¥3, tases V2 ee on Oo +42, 
Vv t 4 t i a F) 
n + 
: | : 3 
v2 4 eo tb2q°v3,, O35 i er be ae 
a ear +2°$2,°v2, + $2, "v3, tal, tele $3, ‘v3, +43, +3 weg +$2)° V3, | $3, +92, ! 
7 : 2° ° ; 2, 3 ce 3 2 ae A 1.63 2 2 
Pe Vig 42° 2," v2, 4 2, | v3, tal, trl: O3q. V3q 103, ~ V2q +62, v3, "$3, +92, | 
a : 3 i. 22 2... ae : 2 
vl, +2°$2,°¥2, 4 $2, “v3, tal, trl at $3, 1 °¥3, +93, +13 - | v2, +$2,°v3, ° $3, +92, 
vl v2 v2 v3 2 
n n n n mv2 = mean(v2) 
pat ag ee k?b./ = k3 = ay = xrspdy V = = : 
eta rl Ka, 3 n rl n 3 *0 Yo = yespdg x0 Ingo Yo lato mv3 = mean(v3) 
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vi. E Se oe ea f DeVere: AV pia pen aye A. oe 
: oe cua y, = Vy, exp! - a ape m pet OAV y cay. exp! - a Oi 
i —<—<— t +h l 41 





, 2 
tka, 4,° yrepd, 4) ~ Vy, exp; ae Oy 


| a a : 
r Y ate ai r er eae a . 
it ct i Xa +Vx, - Vx,c exp: -ar Ay | | Pe x VX exp 2a. ; aa 
a a ee a oe ee 
“Vy Teens r Ty +a 4Vy, - Vy,"e aa 
| t +] = . | ‘ t Yr yy xe : 3 i as 
a Vy,-exp: -ar° A, tk2b, 4, °' dat, 4, ~ |e Hk8, a yespd, 4) - Vy," exp! ~ Or by 
vN H , t a : : 7 ; 
“Xa +Vx, - Vxtexpj -a Ay) 7 
‘ : a : - m ea exp oe: 
) eect St ee HP YR a, , 
— wa ; A 7. ] c Wee i 1 
(Vioul, yy) iy HOY vy #42) “v3, tal, | “Xion, ay 1 | Ret f tr emlcer av | 
Mvicu2 4 -. v2 v3, +3, +92 | 7. 7 ; ; ia 
ue t 1 702, 1 PS, 402, | Vio sy | iy, heft - expla ALi Vy, | 
Viou3, 44 | ne: | i=} tw : is 
~ : i $3, : *v3 +q3 | | Vou, 44 | ; 1 
93, 5 -¥3, 403, Lt enh cay VR | 
HO | 
L oS Sa 
Parameters for this run: 
D- meanl ny N2 =49 of =7.071068 a =0.25 u =5 
AOU = ymean(vl) lOU - ,{ mean( Vioul ) AOU =1.087788 10U =1.197159 


Determine the difference between the true target position, and the filtered position: 


Fx, = X- Xdia, Fy, = Y.7 Ydis, 
the true target position, and the contact position and the true target position, and the IOU-predicted position: 
Cx = Ing, - Xdis, Cy, = lat, = Ydis, iy, = Xiou, - Xdia, ly, = Yiou, = Ydia, 
The root-mean-square error is found for the data: 
‘2 


Sey: ees 32 : 2 ; 2 ey a 
ener | eer cl ener Sra Creme Mt pee | | 


— 


SFx = 2F2x SFy = 2F2y  SCx = 2C2x SCy = 2F2y SIkx = . 2x Sly = <By 


: | (SFx) , (SFy) |. . f (SCx) | (SCy) }, _ | (Silx) | (Sly) }, = = 
rer = | iam (Fx) ‘lam(Fy) 5 Cer = ae aS 5S lerr ileal iz) ek) | § Ferr =0.952507 lerr =1.01829 


Cerr = 1.740496 
The mean filtered and contact velocities is found: 





Bilir yy. <2 -2 | | = = 
FV, = y OLY ae CV, = Uf xrepd, |? + yrapd, mean(FV ) =8.574437 mean(CV ) =10.102023 
mFV = mean(FV) 
sdev( FV) =0.773034 sidev(CV) =1.641369 mCV = mean(CV) 


The RMS error in filtered position 1s compared to the RMS error in the 


contact position: 
TR = Cerr- Ferr TR =0.787989 


The difference from the computed error 1s then computed for 


the current mun: 
COMPerr = AOQU- Ferr COMPerr =0.135281 
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Filtered track using the actual algorithm, with mean glimpse interval 0.05 hour and Sensor AOU = 5nm: 








P Al =0.896146 
A2 =1.251034 
A3 =3.39483! 
, A4 =4.572997 
AS =6.685208 
A6 =16.686341 
A7 =31.313094 
8 AS =31.611555 
Ydis, AO =34.12239 
lay -10 
Y, 
7 N2 =49 
-15 
-20 
=25 “3 6 -4 = 0 2 4 6 8 
Xdis, Ing, .X, 
_ XeSty = 0 yealy = 0 
xest, = eens 2+in(rnd(1 )))-cos 2-4-mnd(1): +Xdis, yest, = vi i++(-2°In(mnd(1 )))+ sun 2° + end(! )! +Y dis, 


“Ve, = Vyrexp -a A, | +3, ri xrspd, 44 “ Vxyrexp! ar Ay st 


Vdiffx, = Ve, - Vx, mean( Vdiffx ) =0.133526 


Vdiffy, = yVe,- Vy, 


Filtered track using estimation technique (using full noise values): 
§ 


= Wytemp 0B, 3, 4g yrmmd, ay -Vyyc exp Ay 


mean( Vdiffy ) =0.030425 


Ps 





Estx, = xest, ~ Xdia, 
Esty, = yest, - Ydis, 
0 pa “pe = ° 
Ex, =: Estx, ” Ey, =. Esty, 
SEx = .Ex SEy = 2Ey 
“5 
_ | (SEx) ,(SEy) . 
Ydis, Eee , last(Ex) last(Ey) | ~ 
a Eerr =! .028216 
on =10 Ferr =0.952507 
st 
tee ymean(vi ) =1.087788 
=ES 
-20 
ge ae 





2 4 6 8 


0 
Xdia, »ing, » XCst, 
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Filtered track using the actual algorithm, with mean glimpse interval 0.1 hours and Sensor AOU = 10nm: 





80 Al =2.089507 
A2 =10.470267 
A3 =14.93814 
A4 =17.572047 
60 AS =20.489058 
A6 =22.232939 
A7 =25.150709 
A8& =25.233842 
Ydis, 40 AD =45.284693 
lat, 
v. PA 
a ee N2 =74 
0 
I ETG =10 =5 0 5 10 15 20 
Xdis, ,Ing, ,X, 
Xealy = 0 yest) = 0 
xest, = ela: ,(-2°In(end(1 )))+ cog 2° eemnd(1 ) | +Xdis, yest, = y Vig’ V2: In( end (1 )))* sin, 2° e+ mnd(1 ) | +Y dis, 


Vdiffx, = xVe,- Vx mean( Vdiffx ) =0.056692 Vdiffy, = yVe,~- Vy, mean( Vdiffy ) =0.013961 


Filtered track using estimation technique (using full noise values): 
80 


: 


xest, ~ Xdia, 
Esty, = yest, ~ Ydia, 


60 


_ || (SEx) ,(SEy) .. 
Yin, : Say Il last(Ex ) last( Ey) ; 2 
as Eerr = 1.925021 
lat, 
sae Ferr =2.165015 
yest Ler re 
aril sh ajmean(vl ) =2.072638 





20 
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Filtered track using the actual algorithm, with mean glimpse interval 0.5 hours and Sensor AOU = 10nm: 


Al =9.779778 
A2 =19.746842 
A3 =23.257026 
A4 =23.336796 
A5 =25.985588 
A6 =26.130136 





A7 =27.402486 
A& =28.369044 
AS =40.739811 
N2 =91 
-40 
0 50 100 150 200 250 300 
Xdis, .Ing, »X, 
Xen = 0 yest) = 0 


vie 


as Nn 


>] ae | 
wba 4 (-2° in( md(1 ))) sin} 2° x°md(1 ): +Ydis, 


xest = 4 


*+,(-2+In(md(1)))> cos: 2+ a+ md(1) | +Xdis,, yest, = 

= t) - e Att - e r = Cd a =F . oe e : ‘if - e ; - e ya 

xVe, = Vx, exp aA HS, 44° “pd, 4, Vx, exp! a Oy yVe, = Vy, exp’ a Oy 43, 41°) yespd, 4, Vy, exp! aby 
Vdiffx, = xVe- Vx mean( Vdiffx ) =-0.104991 Vdiffy, = yVe,~- Vy, mean( Vdiffy ) =0.040473 


Filtered track using estimation technique (using full noise values): 


60 
Estx, = xest, ~ Xdia, 
Esty, > yest, - Ydis, 
E = Estx, Ey =! Esty. ' 
40 Oe tC ety 
SEx = 2Ex SEy = LEy 
_ 1 (SEx) | (SEy) :, 
dis, 2° err = iast(Ex) “laa(Ey). 
ig Eerr =3.402076 
a Ferr =2.944767 
yest, mean( vl ) =3.398984 
—— 0 
-20 
-40 





~50 300 
Xdis, . ing, »Xest, 
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Filtered track using the actual algorithm, with mean glimpse interval 0.5 hours and Sensor AOU = 30nm: 


50 Al =3.9155t2 
A2 =6.91616 
A3 =1t8.287011 
A4 =20.200269 
AS =30.204378 
A6 =34.696944 
A7 =35.833538 
A8 =38.376747 
AD =38.669162 





Ydia, 
Bie 50 
Y, 
Re N2 =96 

~100 

he “100 “50 0 50 100 150 

Xdis, ,lng, .X, 
EB) = 0 yest, = 0 
— ———_ : ; — ——_ _. ; 
xest, = Vtg v2: In(end( )))* cos, 2° w-end(1) | +Xdis, yest, = (Vig? 12° In(end (1 )))- sin, 2-aemd(1)  +¥dis, 


xVe, = Vayrexp a Ay 43, 4y0° xrspd, 44 - Vigrexp! ar Ay!) yVeq = Vy,cexp) -or Ay! He, ao! yrepd, 4, Vyycexp ar Ay 
Vdiffx, = xVe,- Vx, mean( Vdiffx ) =-0.056681 Vdiffy, = yVe,- Vy, mean( Vdiffy ) =-0.145247 


Filtered track using estimation technique (using full noise values): 
50 


Estx, = xest, - Xdis, 

Esty, = yest, - Ydis, 

Ex =: Es a _- 12 
m = Estx, ” Ey, =. Esty, 
SEx = _Ex SEy = cEy 


| (SEx) ,_ (SEy ) We 





Ydis, ESS  last(Ex) last(Ey)! ° 
— Eecr =7.774309 
Ih -50 Ferr =5.817649 
yeat, ymean(vl ) =7.596565 
+_ 

-100 

-150 


-150 -100 0 50 100 150 


0 
Xdia, ,ing, »Xeat, 
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Filtered track using the actual algorithm, with mean glimpse interval | hour and Sensor AOU = 10nm: 


Al =S.02e2727 
A2 =12.437775 
Q3 =17.876472 
A4 =28.662312 
AS =42.644194 
A6 =47.842995 


300 


250 





200 A7 =49.239743 
A8 =50.274019 
AD =53.491883 
Yd 
=" 150 
lat, 
Y, 100 
Bo N2 =45 
50 
0 
-50 
-300 -250 -200 -150 -100 -50 0 50 
Xdis, ,ing, .X, 
xeaty) = 0 yesly = 0 
xest, = Vg i (2° in( end 1)))-cos 2° 4*mnd(! ) +Xdis, yest, = Vig 1(-2° In( end (1! )))> sin, 2° 4+ rnd(! ) | +¥dis, 


xVe, = Weyrexp: ar Ay 43, 4] red, yy Vy exp) arr Ay” Vey = Vy exp! oar Ay | 3, 44°! yrmpdy 44 - Vy enpl -ar° dy 


Vdiffx, = xVe,- Vx, mean( Vdiffx ) =0.472115 Vdiffy, = yVe,- Vy, mean( Vdiffy ) =-0.66914 


Filtered track using estimation technique (using full noise values): 


300 
Estx, = xest, - Xdis, 
Esty, = yest, - Ydis, 
250 

_ 2 _f 2 
Ex, = Estx ” Ey, =. Esty, | 


SEx = Ex SEy = _Ey 
200 


_ PSE) (SE). 


Eerr = |) (2s 
Ydis, 150 TC last(Ex) last(Ey) | 
ins Eerr = 4.03991 

lat, _ 

ae Ferr =3.395647 

yest, 100 ymean(v! ) =3.978631 
= 


50 





~°™~ 300 -250 -200 ~150 ‘100 -50 0 50 
Xdis, -ing, » Xeat, 
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Filtered track using the actual algonthm, with mean glimpse interval | hour and Sensor AOU = 30nm: 


200 Al =0.158067 
A2 =4.15131 
A3 =9.030112 
A4 =14.972823 

lake AS =19212 


A6 =21.285281 
AT =24.284637 





100 A8& =33.048225 
Ydis, A9 =36.179972 
lat, $0 
Y, 
ae N2 =44 
0 
-$0 
50 0 50 100 150 200 
Xdis, Ing, X, 
Xesly = 0 yest) = 0 
xest, = vi, ,(-2+In(rnd(1)))*cos 2-aemnd(1) +Xdis_ yest, = «VE gt vG 2° In(mnd(t )))* sin. 2*a+md(1) +¥dis, 


xVe, = Be exp ca A) +3, ny xrspd, 44 ~Vxyrexp. arr Ay yVe, = Vy, "exp! ar A +3, yrspd, 4, Vy, exp “a Ay 


Vdiffx, = xVe,- Vx mean( Vdiffx ) =-0.495172 Vdiffy, = yVe, - Vy, mean( Vdiffy ) =0.01928 


Filtered track using estimation technique (using full noise values): 


200 
= xest, ~ Xdis, 
= yest, ~ Ydis, 
150 YS’ 2 - 2 
‘ ~ . Ex, = Estx, Ey, =e Esty, 


= Ex SEy = —Ey 


100 
( (SEx) | (SEy) © 
= ‘ .. 
Ydis, , last(Ex) last( Ey), ; 
oa Eerr =9.81617 
a 2 _ Ferr =12.186325 
heal . mean(vl ) =9.549774 
ae 
0 
750 
-100 





Filtered track using the actual algorithm, with mean glimpse interval | hour and Sensor AOU = 60nm: 


Al =0.003175 
A2 =3.366652 
43 =11.418094 
A4 =14.422879 
A5 =19.245857 
A6 =24.664634 
A7 =26.370376 
A8 =30.333605 
A9 =30.916577 


150 





Ydis, 50 
lat, 
Y, 
~~ 9 = 
N2 =54 
-$0 
-100 
100 50 0 50_ 100 150 200 250 
Xdis, sng, .X, 
xen = 0 yest) = 0 
= ———___.. ; — et é 
xest = vin ,(-2+In(md(1 )))* cos: 2° e+ end(1 ) ; +Xdis, yest, = VE VG 2+ In( end (1 )))* sin! 2+ -end(1)  +Y¥dis, 


xVe, = Vxrexp ar dy +3, a xrspd, 44 ~ Vx,sexp' -ar A . pve a Vy," exp! or Dy at +7. Pd, che eee — 


Vdiffx, = xVe,- Vx mean( Vdiffx ) =-0.202094 Vdiffy, = yVe,- Vy, mean( Vdiffy ) =0.294623 


Filtered track using estimation technique (using full noise values): 


150 
Esix, = xest, - Xdia, 


Esty, = yest, ~ Ydis, 





‘Ee? Ee 
100 a Ey, = Ey, 
SEx = _Ex SEy = —Ey 
_ / (SEx) _(SEy ) 2 
Ydis, =e TV, last(Ex) last(Ey), . 
bois Eerr =15.785448 
lat, ” Ferr =12.430599 
yest, ymean(vl ) =15.165514 
+ 0 
-50 
“100 ra “50 0 50 100 150 200 250 
Xdia, , Ing, , xeat, 
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APPENDIX H: KALMAN GAIN VERSUS SENSOR AOU SIZE 


This appendix shows the relationship between the steady state Kalman gain Ka and the Sensor AOU size. 
The impact of using the the approximations from Appendix A is also illustrated. The data for each case was 
written to a seperate file and the graph alone displayed 





2 
ae Oan25 A =:.!l u = 1000 a = .25 B80 4.2 = 25+y2 
2°a 
vleg = rl viey =O v3eq = 3 vig = rl v2y = 0 v39 = 3 
o2 - fe l - exp, -a A o3 -: exp, -a° Aj 
re 4 
= ; ag ; }2 
qi =—: A> =" 2° 1 exp -ar- Ay, - Sl - exp, -@e2* Aisi |! 2 
a os ae q2 = So7f 5 exp[-a Aj + Se expl -a-2° Aj] 
a 
0° 
q3 = Se—— el 1 - exp -a°2°A | 
a 
+2: 62- ’ 27 lr aa? |r 2: °O3 aie 
re a 2°GL-v2, + b2 | ‘v3, tal 7 63, v3, +93 +3 j-: | v2, +6 V3q 9) +92 | ! 
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Using the transition and process noise approximations from Appendix A to simplify the calculations: 
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Kalman gain versus sensor AOU in nm for an interarrival me of 0.01 hours: 
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Kalman gain versus sensor AOU in nm for an interarrival time of 0.1 hours: 
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Kalman gain versus sensor AOU in nm for an interarrival time of | hour: 





KAOU, 


Kalman gain versus sensor AOU in nm for an interarrival time of 2 hours: 
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Kalman gain versus sensor AOU in nm for an interarrival time of 4 hours: 
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